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ABSTRACT 



Providing a simplified model of real terrain has applications to route planning for 
robotic vehicles and militar>' maneuvers. In this thesis I explore planar-patch surface 
modeling to represent terrain in a simple and effective way. In planar-patch surface 
modeling the terrain is subdivided into a set of planar subregions. The homogeneity of 
the gradient within a planar subregion simplifies calculating the cost of traversing the 
region, thus simplifying route planning. I have explored three main strategies to model 
the surface: joint top-down and and bottom-up, strict bottom-up, and presmootliing 
bottom-up approaches. Results of the algorithms are shown graphically by using the 
APL and Grafstat packages, verifying their correctness and accuracy. 
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I. INTRODUCTION 



A. BACKGROUND 

There are many ways to represent a three-dimensional object, usually involving 
complicated mathematical and graphical theories. Representing and displaying a three- 
dimensional object like a terrain is important in applications such as cartography, route 
planning for robotic vehicles, road construction, and planning of military maneuvers. 
A simplified model of real terrain helps the user better understand and analyze it. 

The simplest approximations are piecewise planar, and the very simplest is a 
polyhedron with triangular faces. Assume that a terrain is approximated by a two- 
dimensional array of evenly spaced elevations. Each group of three adjacent points in 
the array defines a planar polyhedral face, and a set of triplets constitutes a planar-patch 
approximation of the surface. Unfortunately, this approach creates many small triangu- 
lar regions, which can make the searching process of route planning too complicated. 

In this thesis, 1 developed three algorithms for planar-patch modeling of evenly 
spaced elevation data, and each has its own advantages and disadvantages. I have used 
the least-squares method to fit a plane over the points of a terrain, a quadtree algorithm 
to subdivide regions, a gradient clustering method to build regions, and smoothing al- 
gorithms on elevation data to improve continuity of adjoining regions. 

One major advantage of planar-patch modeling over other three-dimensional surface 
modeling is that terrain within a patch is homogeneous in the sense of a constant gra- 
dient. This homogeneity simplifies calculating the cost or speed of traversing the region, 
thus making it easier to analyze for route planning. 

Several thresholds affect my algorithms, and different threshold values will result in 
different models. These include the standard deviation of a fitted plane from given 
points, the ratio of magnitude between adjacent points, and the root mean square dif- 
ference of coefficients of planar equations between adjacent regions. The final result of 
this thesis will be provided to Ron Ross (Ph.D. student. Computer Science Department, 
Naval Postgraduate School) for his research in route planning for robotic vehicles. 

B. ORGANIZATION 

Chapter 2 introduces terrain modeling techniques, including least-squares planar 
fitting and cell classification. 
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Chapter 3 is the main chapter of this thesis, defining the planar-patch terrain model. 
Joint top-down and bottom- up and strict bottom-up planar-patch terrain modeling 
strategies are described in detail. Programming techniques such as the quadtree method 
used in the joint top-down and bottom-up algorithm and the region-growing method 
used in the strict bottom-up algorithm are explained. Finally, the issue of continuity of 
adjacent planar-patches along their common boundaries is discussed. Attempts to solve 
this problem, such as smoothing of original elevation data and grouping of the same 
type of points into one region are explained. 

In Chapter 4 comparisons are made between the tw’o algorithms in terms of accu- 
racy and simplicity of terrain representation. Special features and the limitations inher- 
ent to each algorithm are discussed. Experimental results of joint top-down and 
bottom-up, strict bottom-up, and presmoothing bottom-up algorithms are presented in 
graphical form in this chapter. Pictures of the modeled terrain are shown in three di- 
mensions using pictures drawn by the Grafstat and APL packages. Pictures of the 
boundaries of subregions are also included. 

Finally, Chapter 5 summarizes the contributions made by this research and discusses 
some of the possible areas for further research based on this work. The Pascal source 
program is in the Appendix. 
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II. REVIEW OF TERRAIN MODELING 



The most common method of terrain representation is by a set of functional ex- 
pressions such as 

y(arj.') = x^ + 

where f(x,y) is the elevation function for terrain and 
X and y are the latitude and longitude coordinates. 

This permits the simple calculation of the elevation given any combination of x and y. 
It is known that any continuous surface can be approximated with an arbitrarily small 
error by a polynomial of sufficiently high degree. Such polynomials can be derived by 
least-squares surface-fitting methods. 

Surfaces may be defined in tabular form where the value of f is given for a selection 
of representative arguments. To find an arbitrary value of f, a table lookup can be per- 
formed to determine the value of the nearest known points and use them to approximate 
the desired value by interpolation. [Ref 1; p. 212] 

It is possible to specify a surface in terms of guiding points by a generalization of 
the Bezier polynomials or the B-splines. The Bezier polynomial or B-spline method finds 
an approximate curve that passes near a set of given points. [Ref 2: p. 247] 

A. TERRAIN CELL CLASSIFICATION 

I. QUADRATIC APPROXIMATION OF SURFACE BY LEAST SQUARES 

FIT 

Fitting a least-squares quadratic to a 3 by 3 square of points on a surface is done in the 
following way. The fit as a function of x and y is: [Ref 3: p. 55] 

J{x^) = + k2X + + k^x^ 4 - k^xy + k(y>^ 

The square of the error between the fit equation and the given terrain points is: 

+ k^x^k^xy + k^y^ ~Ax,y)f 

The partial derivatives of the squared error with respect to each constant are: 
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To minimize the squared error, each of the above equations is set to zero and simplified. 
Those equations can be further simplified by assuming that the coordinates of each of 
the points in the 3 by 3 matrix, except the center point, is + one distance unit from the 
center point. By the above assumption, any terms in above equations which contain a 
summation of an x or y coordinate with an odd power will go to zero. Simplified 
equations are; 
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0 = ^5 X 4 - 



(15) 




(16) 



From above equations, if we solve for k^, k^, Aj, k^, Aj and A«, we get a quadratic 
equation fitted to the 3 by 3 square of points. 

2. TERRAIN CELL CLASSIFICATION BASED ON QUADRATIC 
APPROXIMATIONS 

The eigenvalues of the Hessian matrix for the quadratic function approximated 
for the 3 by 3 points on the surface are: [Ref 3: p. 61] 



Using the signs of the eigenvalues, we can classify the central point by Table 1 
on page 6. [Ref 3: p. 69] 



{2k, + 2k,) ± v''( -2k, - 2k,f - A{2k,2k, - kl) 
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Table 1. CLASSIFICATION OF SURFACE 
POINTS BY EIGENVALUES OF THE 
HESSIAN MATRIX 



sign of Z, 


sign of 2-2 




- 


0 


+ 


+ 


saddle 


impossible 


depression 


0 


ridge 


plane 


valley 


- 


hill 


impossible 


pass 



B. PLANE APPROXIMATION OF A SURFACE BY LEAST SQUARES FIT 

The elevation z of any point (x,y) can be expressed as: 

The plane that we want to fit over the surface is: 
z = ax ->r by + c 

The square of the error between the fitted plane and the real terrain is: 



y 

The partial derivatives of the squared error with respect to a,b and c are: 



dE^ 

da 
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To minimize the squared error, each of the above equations is expanded and then the 
partials are set to zero. 



0 = “ YYYJ{x^v)x 



X y 



X y 



X y 



X y 



(24) 



0 = Y^axy + 



X y 



X y ^ y ^ y 



(25) 




X y y y ^ y 



(26) 



Here, we have three unknowns, a,b and c, and three equations. So we can solve for 
them. 

C. THE GRADIENT AT A POINT 

One of the primary' criterion we will use for grouping points into a region is the 
gradient. The gradient of a point is a vector that has a magnitude (slope) and a direction 
for its components. We can represent a three by three elevation matrix of terrain points 
as in Table 2 on page 8, where x and y are the coordinates of a point on a terrain and 
f(x,y) is the elevation of the point. 
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Table 2 . A MATRIX OF TERRAIN 



POINTS 



fi;x-l,y+ 1) 


Rx,y+1) 


f{x+ l,y+ 1) 


f(x-l,y) 


f(x,y) 


flx+lj) 


flx-l,y-l) 


Rx.y-l) 


f(x + 1 ,y-l ) 



The gradient vector of a point (x.y) can be approximated in three ways. The first 
method is: 

^i=A-^+ ij') -/w) 



^2 =Ax^+ 1 ) ~Ax,y) 



magniiude{x^v) = \/ [A] + A^) 



-1 ^2 

direciion{x^) = tan ( ) 



The second method is: 

^1 =/('^+ IvK) 

^^2 =Ax^+ 1) -Axy- 1) 

The magnitude(x,y) and direction(x,y) are the same as above. 
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The third way of approximating gradient is by first fitting a quadratic to the three 
by three square of points and then using the coefficients of the quadratic : [Ref 3: p. 60] 

magni[ude{.x^-) = y' (^2 + ^' 3 ) 



direciion{x^) = tan '( ) 



where k-^ and /Cj are the coefficients of the quadratic described in section A. 

The first method considers only one quadrant of the 3 by 3 matrix, thus is biased. 
The second method uses a point from each of the quadrants, and is better than the first 
method in the sense that it utilizes more unbiased information. However, it gives an 
incorrect value in the case as shown in Table 3. 



Table 3. AN EXAMPLE OF A TER- 
RAIN POINTS MATRIX 



5 


2 


3 


J 


1 


3 


-> 


2 


4 



When we calculate the gradient of center point using the second method, it will give 
the magnitude value zero and direction value 45 degree, which is wrong. The third 
method does not have such problem, since it utilizes all eight neighboring points. In this 
case it gives the magnitude value 0.2357 and direction value 45 degree. I used the third 
method in my program. 



9 



The gradient vector of a terrain point points in the x-y direction of maximum slope 
of that point; magnitude ^x,y) is the maximum slope value and direciion(x,y) is the di- 
rection of that slope in the x-y plane. 
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in. PLANAR-PATCH TERRAIN MODELS 



A planar-patch terrain model is a model that represents terrain in the form of a 
polyhedron. This model is different from the common planar terrain model used in 
graphics in that size of a planar-patch is not fixed, and the shape of a patch does not 
have to be regular. 

A. IMPLEMENTATION 

I chose Pascal as an implementation language for this research because of my fa- 
miliarity and experience, and the Waterloo Pascal on the IBM 370/3033AP mainframe 
was the Pascal compiler that was used. 

The APL language and the Grafstat graphics software package were used to verify 
the correctness and accuracy of the algorithms. Even though 1 chose Pascal as my im- 
plementation language, the matrix form of the elevation data suggests APL. With the 
Grafstat software package which runs in the APL environment, one can check the 
intermediate results of an APL program by displaying pictures. Since I experimented 
with different threshold values in my algorithm, APL and Grafstat were useful. The 
three-dimensional pictures in chapter 4 were done with Grafstat. 

Two arrays of records are used as main data structures by all my algorithms: 
point_array and subregion_array. The point_array is a two-dimensional array which 
corresponds to the points in the real terrain. Each element of the array is a record which 
has the elevation, the subregion ID, and the magnitude of the terrain point as its fields. 
From the values of elevation and magnitude, the subregion ID is determined and as- 
signed to the point. The subregion_array is a one-dimensional array of records which is 
created during the subregion creation phase. The indices of this array are the ID's of 
subregions. It keeps necessary information about the subregions such as the equation 
of the plane fitted to the subregion, the boundary of the subregion, and other subsidian.' 
information. 

The general structure of my first two modeling algorithms is shown in Figure 1 on 
page 12. 
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Figure 1. Block Diagram of the First Two Program Approaches 



B. JOINT TOP-DOWN AND BOTTOM-UP PLANAR-PATCH TERRAIN 
MODELING 

My first algorithm consists of two parts. The first part is the top-down subdivision 
of the original elevation data matrix into subregions and the second part is the 
bottom-up merging of adjoining subregions when adjoining subregions have similar 
plane equations. 
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1. TOP-DOWN PHASE 



A quadtree region subdivision method is used during the top-down phase of the 
firtst algorithm. The method generates a quadtree by recursively dividing a two- 
dimensional region into quadrants. Each node in the quadtree has four data elements, 
one for each of the quadrants in the region. [Ref 4: p. 215] The original region is given 
some standard tests such as comparison to a threshold of the maximum magnitude or 
the minimum magnitude of the points within it, depending on the application. Only if 
a region fails the tests and is not of minimum size, it is further subdivided into four su- 
bregions. If a region has both length and width four cells or more, it is not of the min- 
imum size. See Figure 4 on page 16 for the different kinds of minimal subregions. The 
subdivision process described above is repeated for each subregion. Figure 2 on page 
14 shows graphically this algorithm. 

The code works whether a region has a dimensions of odd or even numbers. For 
example, consider 4 by 4 and 5 by 5 matrix regions to subdivide; Figure 4 on page 16 
shotvs how this method works. 

2. THRESHOLDS USED 

Three thresholds are used in the top-down phase: low_bound, upper_bound and 
standard_deviation (SD). 

The threshold upper_bound is set by the maximum climbing capability of the 
vehicle in cases where the polyhedral model will be used for route planning. If the ab- 
solute value of the minimum slope of a region is greater than the upper_bound, which 
means every cell has a slope greater than the upper_bound, then the region is assumed 
being untraversable by the vehicle, and the region need not be subdivided further be- 
cause a good polyhedral fit just doesn't matter for that region. The threshold low_bound 
means that when the absolute value of the maximum slope of a region is less than the 
low_bound, which means every cell has a slope less than the lo'w_bound, then the slopes 
in the region are unimportant. Such regions need not be subdivided either. 

For the region which has maximum slope greater than the low_bound and mini- 
mum slope less than the upper_boiind, the region is subdivided in four parts according 
to the quadtree algorithm to better isolate the unlevel slopes. Then we fit a plane to the 
region and calculate the standard deviation of the plane from the real terrain points. 

The standard deviation (SD) is; 
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QUADTREE REGION SUBDIVISION METHOD 




AN EXAMPLE OF SUBREGIONS CREATED BY THE QUADTREE METHOD 



Figure 2. Quadtree Region Subdivision Method 
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PROCEDURE MAKE_SUBREGIONS(YL,YU,XL,XU : INTEGER; 

prev_region : region_ptr); 

(X YL : LOWER BOUND IN Y AXIS, YU : UPPER BOUND IN Y AXIS, 

XL : LOWER BOUND IN X AXIS, XU : UPPER BOUND IN X AXIS X) 

VAR 

P : REGION^PTR; 

BEGIN 

IF ((YU-YL) > 2) AND ((XU-XL) > 2) THEN (X NOT MINIMUM SIZE x) 

if ( p3 . maxima g_of_regi on > low_bound) then (x test X) 
begin 

nake^subregi ons( C ( yu-yl ) div 2)+yl+l , yu, 

(Cxu-xl) div 2)+xl+l , XU , p); (x 1st quadrant X) 
make_subregions( yl , ((yu-yl) div 2)+yl, 

((xu-xl) div 2)+xl+l , XU , p); (X 2nd quadrant X) 
make_subregions(yl , ((yu-yl) div 2)+yl, 

xl , ((xu-xl) div 2)+xl, p); (X 3rd quadrant x) 

make_subregions( ( (yu-yl ) div 2)+yl+l , yu, 

xl , ((xu-xl) div 2)+xl , p); (x <^th quadrant x) 

P := prev_region ; (x goto parent region x) 

end 

else stop_subdi vi de (X region passed standard test X) 
else stop_subdi vi de (X region has minimum size X) 

END 



Figure 3. Pascal Code for the Quadtree Region Subdivision Method 



where is the elevation of point (x,y) 

ax by -T c is-the plane equation of the region 
n is the total number of points in that region. 

The region is subdivided only if the calculated SD is bigger than the threshold 
SD and it does not have the minimum size. 

3. BOTTOM-UP PHASE 

We now merge adjoining regions having similar plane expressions. The reason 
is that in the process of applying the quadtree subdivision, we often subdivide a region 
which should not have been so completely subdivided. For example, say a region has 
ver\’ steep slopes in the first quadrant and the rest of the quadrants are flat. Since the 
original region will fail in the maximum slope test due to the slopes in the first quadrant, 
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DIFFERENT KINDS OF MINIMUM SIZE 
SUBREGIONS 
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Figure 4. Minitnum-Size Subregions Created by the Quadtree Method; This pic- 
ture shows quadtree subdivision method applied to the regions which 
have dimensions of odd or even numbers, and also shows boundaries 
and different kinds of minimum size subregions. 



the top-down phase will subdivide it in four. However, the first, second and third 
quadrants are all flat and should not be separate. 

The adjoining regions of a region are found from the point_array by looking up 
the subregion ID field of adjoining point. The bottom-up phase is done in two sub- 
phases: row_conibine and col_combine. The rou_combine finds and merges adjoining re- 
gions which are located in a same row, and the col_coninine finds and merges the regions 
in a same column. 
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Figure 5. Subregions after the Merging Process: liicse subregions are created 

by performing the merging process over the initial subregions created by 
the quadtree algorithm. 

To merge adjoining subregions which have similar plane equations, a value 
equation_difference is calculated from the two plane equations and compared to the 
threshold abc_difference. The equation_difference is: 

equal ion jdifference = (Oj — ^2)’ (^1 ~ ^2) + (^1 “ ^2)* 

where a,, d, and c, are the coefficients of plane^ and 
m. and c, are the coefficients of plane^. 

If this is less than the threshold then the two adjoining regions pass the test. 
Then we fit a plane over the two subregions and calculate the standard deviation of the 
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plane, and only if it is less than the threshold standard_deviation we do merge the two 
subregions. 

C. STRICT BOTTOM-UP PLANAR-PATCH TERRAIN MODELING 

1. POINT-CLASSIFICATION PHASE 

In strict bottom-up modeling, w’e first classify every point in a region by at- 
taching two tags to it. The first tag classifies slope of the point: level, safe-slope or 
unsafe-slope. A threshold low_bound defines the low limit of slope and another threshold 
upper_bound defines the upper limit. The second tag classifies the curvature of the point, 
by the terrain cell classification criteria described in chapter 2. The designators are hill, 
depression, plane and special. The tag special refers to all the rest of the cell classifica- 
tions except the above three. A threshold called cun'ature is compared to the eigenvalue 
of the Hessian matrix; if its absolute value is less than the threshold, then the eigenvalue 
is regarded as zero. [Ref 3: p. 69] We group the contuguous points with the same first 
and second tag to form subregions. 

Next the thresholds magnitude_ratio and abc_difference are used to combine the 
subregions. The magnitude_ratio, is used in the region-growing phase and the 
abc_difference is used as in joint top-down and mottom-up modeling. 

2. REGION-GROWING PHASE 

The region-growing phase consists of tw'o subphases: one-dimensional and 
two-dimensional. 

a. ONE-DIMENSIONAL REGION-GROWING 

The one-dimensional phase starts with the first point (/?n) in the first row 
of the point_array. The point becomes the head of a linked list and the subregion_ID 
field of the point is given the value one. The subregion ID starts from one and is in- 
cremented every time a new subregion is created. Whenever a new subregion is created, 
the necessary’ information about the subregion, such as the coordinates of the starting 
point and boolean value that says whether the subregion is active or not, is stored in the 
subregion_array using the subregion ID as an index. Then it examines the next point ( 

/7p) in the same row. If the ratio computed using magnitudes of gradient ( 
aos{magnhudei ,(2 — magnitude 

^ — \ — ) is less than the threshold magnitude ratio, then they 

magnitude ~ 

are linked together to form a linear subregion (Rn)- That is, the next_point field of the 
first point in the point_array is given the x and y coordinate of the second point and the 
subregion_ID field of the second point is given the same value as the subregion ID of the 
first point. The average magnitude of the gradient of the points within the subregion is 
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Figure 6. One-Dimeiisional Region-Growing 

calculated and is stored in the avg_inag field of the subregion record in the 
subregion_array. Then repeatedly for each next point (/>,,) in the same row. the ratio 
between the average magnitude of the gradient of the subregion and the magnitude of 
the gradient of the new point is compared to niagnitude_ratio. If on this basis new point 
(p,J should not join the subregion, a new region starts from that new point and the new 
point becomes a head of another linked list. This process continues for all rows. See 
Figure 6 for an example of this algorithm. 

b. TWO-DIME^’SIONAL REGION-GROWING 

The two-dimensional subphase starts with the first linear subregion (/?,,) 
which was created by the one-dimensional subpha.se. Actually, it starts from the first 
point in the first column of the point_array, which belongs to the subregion one. It ex- 
amines the point below in the same column, which belongs to another subregion, to 
decide whether the average magnitude of the gradient of the subregions are similar. The 
information, average magnitudes of the gradient of the two subregions, is retrieved from 
the the subregion_array using the subregion ID of each point. If the ratio is less than 
magnitude_ratio then they are linked together. The tail of the first subregion is linked 
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Figure 7. T^^o-Dinlensional Region-Gro'>ing: This algorithm is applied after the 

one-dimensional region-growing algorithm. 



to the head of the second subregion, the tail of the second subregion is left nil, all the 
subregion IDs of the points within the second subregion are changed to the first subre- 
gion ID, and the second subregion in the subregion_array is marked inactive. When 
those two linear subregions are merged, the resulting subregion is not one-dimensional 
any more, except if both subregions were single points. Then repeatedly for each next 
point in the same column, the ratio between the average magnitudes of the gradient of 
the subregions is compared to the niagnitude_ratio. If a new subregion does not join the 
old subregion, the control is moved to the new subregion and the process starts from 
that new subregion. This process continues for all columns. This is illustrated in 
Figure 7. .After the two-dimensional phase, we fit a plane to each subregion using 
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least-squares. Then the threshold abc_difference is used to merge adjoining regions that 
have similar plane expressions. The merging is done in two steps. The first step is to 
scan each row of the point_array looking for adjacent points in the same row which have 
diflerent subregion IDs. For eveiy pair of adjacent subregions, the plane expressions 
are consulted from the subregion_array by using the ID numbers of the subregions as 
the indices. The equation_difference of the two plane equations are calculated and 
compared to abc_difference. If it less, the two subregions are merged, point of the sec- 
ond subregion to the tail point of the first subregion. The second step is identical except 
that scanning for merging is done by column. 

D. AN ALGORITHM IMPROVING THE CONTINUITY OF ADJOINING 
REGIONS 

One common problem with planar-patch terrain modeling is that planar patches for 
two adjacent subregions seldom yield a surface that is continuous along the line seg- 
ments bisecting the points used to create the regions. Since the initial requirements for 
this research rule out any other terrain modeling than planar-patch, we had to come up 
with some idea to ensure continuity of the patches. One simple solution is to model 
terrain entirely with small triangles, each of which exactly fits three adjacent terrain 
points, as described in Chapter 1. The number of triangles necessary is very' large. 

But a minor modification to the two modeling approaches explored in this research 
can give smoothness. One can create supplementary triangular subregions to fill the gap 
between adjoining subregions. The gap between adjoining subregions have two lines, one 
for each side of the gap. The lines are the edges of the adjacent planes. We can define 
two triangles on the gap by grouping the first end point of the first edge to the second 
edge and second point of the second edge to the first edge. See Figure 8 on page 22. 
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TWO TRIANGLES CONNECTING 




Figure 8. Filling the Boundan with Two Triangles 



This is better than the many-triangle model, but is still not optimal, because the number 
of subregions will be increased considerably by the supplementary triangular subregions. 

We have made a dilTerent attempt to solve this problem using smoothing of the el- 
evation array and an expansion to three dimensions of Rolle's theorem. The basic idea 
is that if we lit planar-patches to an everywhere-convex or evers'vvhere-concave region, 
the fitted patches are much more likely to be continuous. This means we should try to 
ensure that all the points within large regions are classified as either planar points, hill 
points or depression points. See Figure 9 on page 23. Smoothing the elevation data 
will help: other kinds of points will tend to be changed into one of the three. 

Unfortunately, because of time, we were not able to solve this rather complex 
problem completely. It is recommended for future research. 




INTERSECTION 




THE INTERSECTION TENDS TO LIE BETWEEN THE TWO REGIONS 
IF WE FIT PLANES OVER THE EVERY-WHERE CONVEX OR 
EVERY-WHERE CONCAVE TERRAIN 



Figure 9. Hon Surface Convexity Helps 



23 



IV. DISCUSSION 



A. ACCURACY ANALYSIS 

The joint top-down and bottom-up and strict bottom-up approaches produce fairly 
good terrain models in terms of accuracy and simplicity. Here accuracy means whether 
the model represented the real terrain without losing anything significant, and simplicity 
means the number of subregions. 

The overall accuracy of a model can be verified by comparing the pictures of the real 
terrain and model. But there is another thing that describes the accuracy of a model, 
standard deviation. We can evaluate the local accuracy in terms of the standard devi- 
ation of the elevation dilTerence for a point in the real terrain and the corresponding 
point in the model. In the top-down approach, each subregion must have a fit standard 
deviation less than the the threshold standard_deviation. Let's assume that the elevation 
difierence between a point on the real terrain and a point on the model has a normal 
probability distribution. If we know the standard deviation of the elevation differences 
of a subregion: 

1. About 68 percent of the points will have elevation difference less than plus or nii- 
nus one standard deviation. 

2. About 98 percent of the points will have elevation difference less than plus or mi- 
nus two standard deviations. 

3. Practically all points (99.73 percent) 

will have elevation difference less than plus or minus three standard deviations. 

So for example, if we limit the model to have the worst case local elevation differ- 
ence, say 10 feet, with 98 percent certainty, we can set the standard deviation threshold 
to 5 feet. Consequently every subregion will have a standard deviation less than 5 feet, 
and 98 percent of the points will have an elevation difference less than 10 feet. 

The threshold standard deviation was not used in the strict bottom-up approach so 
the above argument does not hold. 

B. LIMITATIONS OF THE MODEL 

In the joint top-down and bottom-up approach, the minimum number of points in 
a subregion is four and the shape of initial subregion must be a rectangle. This re- 
striction is due to the quadtree. These limitations prohibit the approach from picking 
up linear (one-cell-wide) terrain features and creating linear subregions. This limitation, 
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however, does not exist in the strict bottom-up approach owing to the linked list of 
points used to represent a region. When a linear region is a straight horizontal or ver- 
tical line, one cannot calculate a plane equation. 

C. EXPERIMENTAL RESULTS 

We now present the output of our terrain modeling program in graphical form. It is 
applied to sample elevation data acquired from the Fort Hunter Liggett terrain data 
base. Shown are the pictures that are produced by our modeling programs using differ- 
ent thresholds. 

As was mentioned in chapter 3, in joint top-down and bottom-up modeling, the 
thresholds were low_bound and upper_bound for gradient magnitude, standard_deviation, 
and abc_difference. In strict bottom-up modeling, the thresholds were low_bound and 
upper_bound of gradient magnitude, curvature, abc_difference, and inagnitude_ratio. 

Three sets of pictures are given sequentially in the following section. The first set 
is produced by the joint top-down and bottom- up approach (referred as the first ap- 
proach in the following pictures), the next by the strict bottom-up approach (referred 
as the second approach), and the last by the bottom-up approach applied to the 
smoothed data (referred as the third approach). 

In the array pictures, each subregion has a unique code. If a number does not have 
a prefix then the number is actually the subregion ID. If it does have a prefix then the 
number has to be converted. See the conversion example below. 

23 = = > subregion ID 23. 

.23 = = > subregion ID 123. 

:23 = = > subregion ID 223. 

-23 = = > subregion ID 323. 

0 = = > either an outermost boundarx' 
or inside area of a region. 
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1. EXPERIMENTAL RESULTS 1 OF THE FIRST APPROACH 
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Figure 10. Sample Real Terrain: This is the real terrain that were used in these 
experiments. 
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Figure II. Subregions Created by the Quadtree Method in the First 
Approach: Thresholds: lo\v_bound = 0.1; upper_bound = 0.6; 

standard_dev:ation = 3. The above picture shows the initial subregions 
created by quadtree algorithm. 
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Figure 12. Subregions Created after the Merging in tlie First Approach: Thresh- 

olds: abc_difierence= 10; siandard_deviation = 3; Some of the initial 
subregions created by quadtree algorithm are merged by merging 
process, thus yielding less subregions. 
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Figure 13. Terrain Model Created by the First Approach: This is the final prod- 

uct of the first approach. 
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Figure 14. Terrain Model Created by the First Approach Showing Boundaries 
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2. EXPERIMENTAL RESULTS 2 Of FIRST APPROACH 
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15. 




Subregions 


Created 


by 


the 




Quadtree 


Method 


in 


the 




First 





Approach: Thresholds: low_bound = 0.1; upper_bound = 0.6; 

siandard_dc\ iation = 7. The above picture shows initial subregions cre- 
ated by the quadtree algorithm using a different threshold (the 
standard_deviation is changed). 
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Figure 16. Subregions Created after the Merging in the First Approach: Thresh- 

olds: abc_difrcrence= 30; standard_deviation= 7; Some of the initial 
subregions created by quadtree algorithm are merged by the merging 
process; The thresholds standard_de%iation and the abc_dilTerence are 
changed. 
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Figure 17. Terrain Model Created by the First Approach: Final product of the 

first approach using a different set of thresholds. 
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Figure 18. Terrain Model Created by the First Approach Showing 
Boundaries: The planar patches are different from the planar patches 
produced from the previous experiment. 
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3. EXPERIMENTAL RESULTS 1 OF THE SECOND APPROACH 
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Figure 19. Subregions Created by the Region Growing Algorithm in the Second 
Approach: Thresholds: low_bound = 0.1; uppcr_bound = 0.6; 

magnitude_ratio = 0.1. This picture is produced after performing the 
one-dimensional region growing process (row region growing). 
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Figure 20. Subregions Created after the Merging in the Second 
Approach: Threshold: abc_difTerence = 10; the two-dimensional re- 

gion growing algorithm is applied over the subregions created by the 
one-dimensional region growing algorithm and then merging is per- 
formed. 
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Figure 21. Terrain Model Created by the Second Approach: Final product of the 
second approach. 
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Figure 22. Terrain Model Created by the Second Approach Showing Boundaries 
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-1. EXPERIMENTAL RESULTS 2 Of THE SECOND APPROACH 
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Figure 23. Subregions Created by the Region Gro.ting Algorithm in the Second 
Approach: Thresholds: lo\v_bound = 0. 1; upper_bound = 0.6; 

magniiude_ratio = 0.2. This picture is produced after performing onc- 
dimcnsional region growing using a different threshold {the 
magnitudc_ratio is changed). 
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Figure 24. Subregions Created after the Merging in the Second 
Approach: Threshold values used: abc_difference = 20. Two dimen- 

sional region growing is applied and then merging is performed using 
different thresholds (abc_difference and magnitude_ration are 
changed). 
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Figure 25. Terrain Model Created by the Second Approach: Final product of the 

second approach. Notice that this picture is different from the picture 
produced by the previous experiment. 
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Figure 26. Terrain Model Created by the Second Appr 9 ach Showing Boundaries 
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5. EXPERIMENTAL RESULTS OF THE THIRD APPROACH 
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Figure 27. Smoothed Sample Terrain: The original sample terrain is smoothed. 

Bottom-up modeling wall be performed on this smoothed terrain. 
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Figure 28. Subregions Created by the Tliird Approach: Thresholds: 
lo\v_bound = 0.1; upper_bound = 0.6; magnilude_raiio = 0.1. 
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Figure 29. Subregions Created after the Merging in the Tliird 
Approach: Threshold: abc_dinercnce= 10. 
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Z AXIS (FEET) 




Figure 30. Terrain Model Created by the Third Approach: This picture looks 

almost the same as the the picture produced by the second approach; 
this approach produced less subregions than the second approach. 
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Figure 31. Terrain Model Created by the Third Approach Showing Boundaries 
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V. CONCLUSIONS 



In this research, we explored three types of planar-patch terrain modeling. The purpose 
of such models is a simple and yet accurate view of terrain. Since our models were de- 
veloped for route planning, we required a minimal number of subregions, since this di- 
rectly affects the cost of the route planning. 

The first algorithm (joint top-down and bottom-up algorithm) takes the longest ex- 
ecution time among the three. The second (strict bottom-up) and the third (bottom-up 
algorithm applied to the smoothed elevation data) algorithms take almost same time, 
though the third algorithm requires the input elevation data to be smoothed first, re- 
quiring some extra preprocessing time. As far as result accuracy (simplicity) is con- 
cerned, the first algorithm is better than the others. The second and third do not use a 
standard deviation threshold, resulting in possible inaccuracy, and tend to produce more 
subregions. Since all of the three algorithms use almost the same data structure, a two- 
dimensional point_array and a one-dimensional subregion_array, they take almost the 
same space. As for the application environments of these algorithms, the first algorithm 
could be used in situations where the accuracv is critical, while the second and third 
algoritms could be preferred in situations where a short execution time is required and 
linear terrain features should be identified. 

Since the algorithms used in this research operate on large arrays and involve heavy 
computations, a multiple-processor computer architecture would yield superior process- 
ing capability since the same operations could be performed for each part of the input 
matrix simurntaneously. As I conclude this research, I regret that I could not solve the 
continuity issue completely. 
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APPENDIX A. PASCAL CODE OF JOINT TOP-DOWN AND 
BOTTOM-UP ALGORITHM 



(*$s350000*) 

( *^V****VcVf***Vf'>V*^V*^VVc*Vr****Vr*Vr***’>V****^V*Vc**Vf****’*Vf****’'>V';V’>V’jVVf'}V'^**’>V’)V'}V**’jV^V’>V** 



Author 


= Yee, Seung Hee. 


Date 


= 3 


June 1988. 


Input 


= 1. 
2. 


Elevation data acquired from Fort. Hunter Liggett 
data base (40 by 40 square). 

file of magnitude of gradient of the above elevation 
data created by using the program in Appendix 3. 


Output 


= 1. 
2. 


Elevation data of the model created by this program. 
Elevation data of the boundary points, other points 
are given the value one so that when we display 
this model in three-dimensional pictures, we can 
only see the boundaries of the patches. 



Computer = Waterloo Pascal on IBM 370/3033AP, 
Naval Postgraduate School. 
Description = 



This program is the Pascal implementation of the 
Joint top-down and bottom-up terrain modeling. 



) 



program top_down( input , output); 
const 

mtinft = 3. 2808; 

low_bound = 0. 1; (* mag > 10 % slope then divide *) 

(* mag <= 10 % then disregard *) 
upper_bound = 0. 6; (* mag < 60 % then do not subdivide *) 

max_sd = 3; 
abc_dif f erence =10 ; 
max_group_num= 300; 
row = 40; 
col = 40; 
increment = 3; 

type 

elrange = 1. . 10000; 
pointrec = record 

el : elrange ; 
mag: real; 

gp : 0. . max_group_num; 
end; 

magrec = record 

mag : real; 
r ,c : 1. . row; 
end; 

point_array=array ( . 1. . row, 1. . col. ) of pointrec; 
grouppointer = -*grouprec ; 
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grouprec = record 

rl.rh, 

cl ,cr : 1. . row; 
magmax,magmin : magrec ; 
up : grouppointer; 
end; 

constarray = array ( . 1. . 3 , 1. . 4. ) of real; 

nodepointer = •’node; 
node = record 

r,c : integer; 
next : nodepointer; 
end; 

coefrec = record 

crl ,crh, 

ccl ,ccr : 1. . row; 

a,b,c : real ; 
minmag,maxmag : real; 
sd : real ; 
ee,ex,ey : real; 
side : 0. . raax_group_num; 
constarr : constarray; 
active : boolean; 
list : nodepointer ; 

end; 

subregion_array = array ( . 1. . max_group_num. ) of coefrec ; 



var 

points : point_array ; 
groups : grouprec ; 
subregion : subregion_array ; 
ok : constarray; 
data,datam,doutl,dout2 : text; 
counter, r,c,i : integer; 
avg,ee,ex,ey : real; 
top,p ; grouppointer; 

temparr: array ( . 1. . max_group_num. ) of boolean; 

function variance (a,b,c,ee,ex,ey: real; ok: constarray): real; 
begin 

variance := ( ee-( 2*a*ex) -( 2*b*ey) -( 2*c*ok( . 3 ,4. ) ) 
+(a*a*ok(. 1,1. ))+(2*a*b*ok(. 2,1. )) 
+(2*b*c*ok(. 2,3. ))+(2*a*c*ok(. 1,3. )) 
+(b*b*ok(. 2,2. ))+(c*c*ok(. 3,3. ))) 
/(ok(.3,3. )-l); 

end; 

procedure initialize; 
var 

r,c : integer; 
begin 

counter := 0 ; 

for c : = 1 to col do 
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for r : = 1 to row do 

points(. r. )(• c. )• gp := 0 ; 

for r : = 1 to max_group_num do 
begin 

subregion( . r. ) • a “999; 
subregion( . r. ) • b := -999; 
subregion(. r. ). c := -999; 
subregionC . r. )• active := false; 
subregion( . r. )• side := 0; 
subregion(. r. )• sd := -999; 

end; 

end; (* procedure initialize *) 



procedure 

var 



getconst(var points: point_array; rl,rh,cl,cr: integer; 
var arr: constarray; var ee,ex,ey :real); 



i, j: integer; 



begin 

for i : = 1 to 3 do 
begin 

for j ;= 1 to 4 do 
arr(. i, j. )= =0; 

end; 



ee : = 


0; 
















ex : = 


0; 
















ey : = 


0; 
















for i 


: = 


rl to 


rh do 










begin 




















for 


j • = 


cl to 


cr do 










beg 


in 


















arr(. 


1,1. 


) 


=arr(. 


1,1. 


)+( i*i); 








arr( . 


1,2. 


) 


=arr( . 


1,2. 


)+(j*i); 








arr( . 


1,3. 


) 


=arr( . 


1,3. 


)+ i; 








arr( . 


1,4. 


) 


=arr( . 


1,4. 


)+(points(. i, j 


. ).el*i); 






arr( . 


2,1. 


) 


=arr( . 


2,1. 


)+U*j); 








arr( . 


2,2. 


) 


=arr( . 


2,2. 


)+(j*j); 








arr( . 


2,3. 


) 


=arr( . 


2,3. 


)+ j; 








arr( . 


2,4. 


) 


=arr(. 


2,4. 


)+(points( . i , j 


. ).el*j); 






arr( . 


3,1. 


) 


=arr(. 


3,1. 


)+i; 








arr(. 


3,2. 


) 


=arr( . 


3,2. 


)+j; 








arr( . 


3,3. 


) 


=arr( . 


3,3. 


)+i; 








arr(. 


3,4. 


) 


=arr( . 


3,4. 


)+points(. i, j. 


).el; 






ee : = 


= ee 


+ 


sqr(points(. i, j. ). el); 








ex : = 


= ex 


+ 


(points( . i 


., j. ). el * i ); 








ey : = 


= ey 


+ 


(points(. i 


., j. ). el * j ); 





end (* j *) 

end; (* i *) 

end; (* procedure get const *) 
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procedure gauss(ok: constarray; var aa,bb,cc; real); 



var 

i.j : integer; 

temp : array (.1..4. ) of real; 

t ; real; (* error checking purpose *) 

begin 

i := 1; 

while ok(.i,l. ) = 0. 0 do i := i + 1 ; (* row rotation *) 

for j ;= 1 to 4 do 

begin 

temp(. j. ) := ok(. 1, j. ); 
ok(. 1, j. ) := ok(. i, j. ); 
ok(. i,j. ) := temp(. j. ); 
end; 

t := ok(. 1,4. ); 

for i : = 2 to 3 do 

for j : = 4 downto 1 do 

ok( . i , j . ) : = ok( . i , j . )*ok( .1,1.) + ok( . 1 , j . )*( -ok( . i , 1. ) ) ; 
i : = 2; 

while ok(.i,2. ) = 0. Odo i:=i + l; 

for j := 1 to 4 do 

begin 

temp(. j. ) := ok(. 2, j. ); 
ok(. 2, j. ) := ok(. i, j. ); 
ok(. i, j. ) ; = temp(. j. ); 
end; 

for j : = 4 downto 1 do 

ok(.3,j.) := ok(.3,j. )*ok(.2,2. ) + ok( . 2 , j. )*( -ok( . 3 , 2. ) ); 

if ok(.3,3. ) <> 0.0 then 
begin 

ok(. 3,4. ) ;= ok(. 3,4. ) / ok(. 3,3. ); 

ok(.2,4. ) ;= ( ok(.2,4. ) - ok( . 3 ,4. )*ok( . 2 , 3. ) ) / ok(.2,2. ); 
ok(.l,4. ) ;= ( ok(.l,4. ) - ok(. 2,4. )*ok(. 1,2. ) 

- ok(.3,4. )*ok(. 1,3. )) / ok(.l,l. ) ; 
aa := ok(.l,4. ); (* solution *) 

bb := ok(.2,4. ); (* solution *) 

cc := ok(.3,4. ); (* solution *) 

(* gauss error = t-( aa*ok( . 1 , 1. )+bb*ok( . 1 , 2. )+ cc*ok( . 1 , 3. ) ) *) 
end; 

end; (* procedure gauss *) 



procedure prefix(num: integer); 
begin 

case ((num div 100) mod 10) of 

0 : write( ' ' ); 

1 : write( ' . ' ); 

2 : write( ' ; ' ); 

3 : write( ' ~ ' ); 

4 : write( '+' ); 
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5 : write( 

6 : write( 

7 : write( 

8 : write( ^ ? ^ ); 

9 : write(^=*); 

end; 

end; (* procedure prefix *) 

procedure printsubl(points; point_array); 
var 

: integer; 

begin 

page; 

for k : = 1 to counter do 
with subregionC . k. ) do 

for i := crl+1 to crh-1 do 
for j := ccl+1 to ccr-1 do 
points(. i, j. )• gP : = 0; 
for r ; = 1 to row do 
for c : = 1 to col do 
begin 

if c<> col then begin 

pref ix(points( . r,c. ). gp); 

write((points(. r,c. )• gP 100): 2); 

end 

else begin 

pref ix( points ( . r ,c. )• gp); 

writeln( (points( . r ,c. )• gp 100): 2) 

writeln; 

end; 

end; 

end; (* procedure printsubl *) 



procedure printsub2(points: point_array); 
var 

k,r,c,rl,r2,cl,c2 : integer; 
current : nodepointer; 
begin 

for r:= 2 to row-1 do 
for c: = 2 to col-1 do 

points(. r,c. )• gp : = 0 ; 
for k: = 1 to counter do 
if subregionC . k. )• active then 

begin 

current := subregion( , k. ) • list ; 
while current <> nil do begin 
rl: =current^. r; 
cl: =current-*. c; 

if current^, next <> nil then begin 
r2: =current^. next-', r; 
c2; =current^. next-*, c; 

end 

else begin 
r2 := rl; 



53 



c2 : = cl; 
end; 

if rl > r2 then begin r: = r2; r2: = rl; rl:= r; end 



if cl > c2 then begin c: = c2; c2: = cl; cl:= c; end 

if rl=r2 then for c : = cl to c2 do 

points(. rl,c. ). gp ;= k; 

if cl=c2 then for r : = rl to r2 do 

points(. r ,cl. )• gp k; 

current := current-’, next; 
end; 
end; 

for r : = 1 to row do 
for c : = 1 to col do 
begin 

if c<> col then begin 

pref ix(points( . r,c. ). gp); 
write((points(. r,c. ). gp mod 100): 2); 
end 

else begin 

pref ix(points( . r,c. ). gp); 
writeln((points(. r,c. ). gp mod 100): 2); 
writeln; 
end; 

end; 

end; 

(* procedure printsub2 *) 



procedure get_nodes( var subregion: subregion__array); 
type 

direction = ( up , down, right , left ); 

var 

i,j,k : integer; 

current 5 head, p : nodepointer; 
temp : direction ; 

procedure get_next(var p : nodepointer; 

var temp : direction); 

var 

i,j : integer; 
done : boolean ; 



procedure do_dir(var done: boolean ; dir: direction ); 
begin 

new(p); 
p-. r := i ; 

P"- c := j ; 
done := true; 
case dir of 

up : temp : = up; 

down : temp : = down; 

right : temp := right; 
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left : temp := left ; 

end; (* case 

end; (* sub_sub_ procedure do_dir *) 



begin (* sub procedure get next *) 
done := false ; 
i ;= P-. r; 
j := c; 
case temp of 

right : begin 

j := j+1; 

while not done do begin 

if points(. i-1, j. )• gP = k 
then do_dir(done,up); 
if (not done and 

(points(. i, j + 1. )• gP k)) 

then do_dir( done , down) ; 
if( (subregion( . k. ). list^. r = i) 

and ( subregion( . k. ). list-*, c = j)) then 
begin 

new(p); 

P-. r : = i ; 

P-. c : = j ; 
done : - true; 
end; 

if not done then j:= j+1; 
end; (* while *) 
end; 

left : begin 

j := j"l; 

while not done do begin 

if points(. i+1, j. ). gp = k 
then do_dir( done, down); 
if (not done and 

(points(. i, j-1. ). gp <> k)) 
then do_dir(done,up); 
if not done then j:= j-1; 
end; (* while *) 
end; 

up : begin 

i := i-1; 

while not done do begin 

if points(. i, j-1. ). gp = k 
then do_dir(done, left); 
if (not done and 

(points(. i-1, j. ). gp <> k)) 
then do_dir( done, right); 
if not done then i: = i-1; 
end; (* while *) 
end; 

down : begin 

i : = i+1; 

while not done do begin 

if point s( .i,j+l. ).gp = k 
then do_dir( done, right); 
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if (not done and 

(points(. i+1, j. )« gP <> k)) 

then do_dir( done, left); 
if not done then i: = i+1; 
end; (* while *) 
end; 

end; (* case *) 

end; (* sub procedure get_next *) 



begin (* procedure get nodes*) 
for k : = 1 to counter do 
if subregionC . k. )• active then 
begin 

new(p); 

p-'. r : =subregion( . k. ) • cr 1 ; 
p-*. c : =subregion( . k. )• ccl ; 
while ( points( . p-'. r-l,p-*. c. )• gP = k) do 

p-*. r := p-*.r-l; (* move to the upper most row of the gp 

while ( points( . p-'. r jp-'. c-1. ). gp = k) do 

p-*. c := p-*. c-1; (* move to the left most col of the gp 

subregion( . k. )• list : = p ; 
current : = p ; 
temp : = right ; 
repeat 

get_nex t ( p , t emp ) ; 
current-*, next : = p ; 
current : = p ; 

until (subregion(. k. ). list--, r = p"*. r) 

and ( subregion( . k. ), list-*, c = P“*. c) ; 
p-*. next : = nil; 
end; 



(* for k : = 1 to counter do 
if subregion( . k. ). active then 
begin 

writeC’gp ’ ,k: 3, ***=>' ); 
current := subregion( . k. ) . list ; 
while current <> nil do begin 

write( current-*, r: 2, ’ , ’ , current-*, c: 2, ’=>’ ); 
current : = current-*, next; 
end; 

writeln; 

end*) 

end; (* procedure get^nodes *) 



procedure makesub( r 1 , rh, cl , cr .-integer; prev: grouppointer); 
type 

why = (nondivisible,magnitude,std, steep); 



var 

dif,r,c : integer; 
flag : why; 
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sd,a,b,cc : real; 

procedure stopdivide( f lag : why ); 
var 

r,c : integer; 
begin 



counter : = counter + 1; 
with subregionC . counter. ) do 
begin 

crl := rl; 
crh : = rh; 
ccl : = cl; 
ccr : = cr; 

maxmag : = p-*. magmax. mag ; 
minmag := p“*. magmin. mag ; 
end; 

for c : = cl to cr do 
for r : = rl to rh do 

points(. r. )(. c. ). gp := counter; 

(* case flag of 

nondivisible : 

write('* non divisible area!*); 
magnitude 

write(** stop divide! by mag.*); 
std : 

write(** stop, it fits well.*); 
steep : 

write(** stop, it is too steep.*); 

end; 

writelnC** area r 1 , rh, cl , cr=* , r 1: 3 , rh: 3 , cl: 3 , cr: 3 , 
* has gpnum=* , counter: 3 ); 

*) 

p := prev ; 

end; (* sub procedure stop divide *) 



begin (* procedure makesub *) 
new(p); 

P“*. magmax. mag := -999; 
p-*. magmin. mag := 999; 
p-». rl : = rl; 
p-». rh : = rh; 

P“*. cl : = cl; 
p“*. cr : = cr; 
p-*. up : = prev ; 

for c : = cl to cr do 
for r : = rl to rh do 
begin 

if points( . r. )( . c. ). mag > p-*. magmax. mag then 
begin 

p-’. magmax. mag := points( . r. )( . c. ). mag; 
p“’. magmax. r : = r; 
p-*. magmax. c := c; 
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end; 

if points( . r, ) ( • c. ) • ^ p-’. magmin. mag then 

begin 

p-». magmin. mag := points (. r. )(• c- )• 
p“*. magmin. r : = r; 
p*’. magmin. c : = c; 

end; 

end; 

(* max, min magnitude for region rl,rh,cl,cr before division *) 

get cons t ( point s, r l,rh,cl,cr, ok, ee, ex, ey); 
gauss(ok, a,b,cc); 

sd := sqrt(variance(a,b,cc,ee,ex,ey,ok)); 

(* sd for gp before division *) 

if ((rh - rl) > 2) and ((cr - cl) > 2) then 
if (p-*. magmax. mag > low_bound) and 
(p-*. magmin. mag < upper_bound) then 
if sd > max_sd then 
begin 

makesub(rl , ((rh-rl) div 2) + rl, 

cl , ((cr-cl) div 2) + cl, p); 

makesub( (( rh-rl) div 2) + rl + 1 , rh, 
cl , ((cr-cl) div 2) + cl , p); 

makesub(rl , ((rh-rl) div 2) + rl, 

((cr-cl) div 2) + cl + 1 , cr , p); 
makesub( (( rh-rl) div 2) + rl + 1 , rh, 

((cr-cl) div 2) + cl + 1 , cr , p); 
p := prev ; 

end 

else stopdivide( std) 
else stopdivide(raagnitude) 
else stopdivide(nondivisible); (* smallest area *) 
end; (* procedure makesub *) 

procedure do_combine( var pt: point_array; 

var subregion: subregion_array; 

var last, new : integer); 

var ii,jj : integer; 

tail, current : 0. . max_group_num; 
begin 

subregion(. last, ). active := true; 
subregion( . new. ). active := false; 
current := subregion( , last, ), side ; 
if current = 0 then 
begin 

subregion( . last. ). side := new; 
current : = new; 

end 

else begin 

while current <> 0 do 
begin 

tail := current ; 

current := subregion( . current. ), side ; 

(* traverse *) 
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end; (* subregion arr *) 

subregion(. tail. )• side := new; 

current := new; (* 

end; 

while current <> 0 do (* reassigning gp # 
begin 

for ii := subregion( . current. ). crl to 
subregion(. current. ). crh do 
for jj := subregion( . current. ). ccl to 
subregion(. current. ). ccr do 
pt(. ii, j j. ). gp : = last; 
current : = subregion(. current. ). side; 
end; (* while *) 
for ii : = 1 to 3 do 
for j j := 1 to 4 do 

subregion(. last. ). constarr(. ii, j j. ) : = 

subregion(. last. ). constarr(. ii, jj. ) + 
subregion(. new. ). constarr(. ii, j j. ); 
subregion( . last. ) . ee := coefs(. last. ). ee + 

subregion( . new. ) . ee ; 
subregion( . last. ). ex := coefs(. last. ). ex + 

subregion( . new. ). ex ; 
subregion( . last. ). ey := coef s( . last. ) . ey + 

subregion( . new. ). ey ; 
with subregion(. last. ) do 
begin 

gauss ( cons tar r , a , b , c) ; 

sd := sqrt(variance(a,b,c,ee,ex,ey,constarr)); 

end; 

end; (* sub procedure do_combine *) 



(* to link new gp *) 
into existing link *) 

to new region *) 



procedure combine(var pt: point_array; 

var subregion: subregion_array); 

var 

tabc_dif f erence, i, j , last , new : integer; 
joinabCjjoinsd, joined : boolean; 

procedure compareabc( var join: boolean); 
var 

temp : real ; 
begin 

temp := sqrt( sqr( subregion( . new. ) . a - 

subregion(. last. ). a) + 
sqr ( subr egion( . new. ) . b - 
subregion(. last. ).b) + 
sqr(subregion(. new. ). c - 
subregion(. last. ). c) ); 
if temp <= tabc_dif ference then join := true 
else join : = false; 

end; (* sub procedure compareabc *) 

procedure checksd( first , second: coefrec; var join: boolean); 
var 
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ii,jj : integer; 
begin 

for ii : = 1 to 3 do 
for j j : = 1 to 4 do 

first. constarr(. ii, j j. ) : = 

first. constarr(. ii, jj. ) + 
second. constarr(. ii, jj. ); 
first.ee := first.ee + second.ee ; 
first, ex := first, ex + second, ex ; 
first, ey := first, ey + second, ey ; 
with first do 
begin 

gauss(constarr ,a,b ,c); 

sd := sqrt(variance(a,b,c,ee,ex,ey,constarr)); 

if sd > max_sd then join := false else join := true ; 
end; 

end; (* sub procedure check_sd_before_combine *) 

begin procedure combine 

tabc_dif ference := increment; 

while tabc_dif ference <= abc_dif ference do begin 

joined : = false; 

for i : = 2 to (row-2) do 

begin 

last : = pt(. i,2. ). gp ; 
for j : = 2 to (col-2) do 
begin 

new : = pt(. i, j + 1. ). gp; 

if ((i=row-2) and (j=col-2)) and (not joined) then 
subregion( . new. ). active := true; 
if (newolast) and 

(subregion(. new. ). maxmag <= upper_bound) and 
( subregion( . last. ) . maxmag <= upper_bound) then 
begin 

compareabc( joinabc) ; 
if joinabc 
then begin 

checksd( subregion( . last. ) , 

subregion(. new. ) , joinsd ); 
if joinsd then begin 

do_combine( pt , 

subregion, last , new); 
joined : = true; 

end 

else begin 

subregion( . last. ). active := true; 
joined := false ; 
last : = new ; 

end; 

end 

else begin 

subregion( . last. ). active := true; 
joined : = false ; 
last : = new ; 

end; 

end 
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else 



begin 

subregion(. last. )• active := true; 
joined := false ; 
last := new ; 

end; 

end; 

end; (* row combine finished 
for j : = 2 to (row-2) do 
begin 

last := pt(. 2, j. ). gp ; 
for i := 2 to (col-2) do 
begin 

new := pt(. i+1, j. ). gp; 
if (newolast) and 

(subregion( . new. )• tninmag <= upper_bound) and 
(subregion(. last. ). rainmag <= upper_bound) then 
begin 

compareabc( joinabc); 
if joinabc 
then begin 

checksd( subregion(. last. ), 

subregion( . new. ),joinsd ); 
if joinsd then begin 

do_combine(pt , 

subregion, last , new); 
joined := true; 

end 

else begin 

subregion( . last. ). active := true; 
joined : = false ; 
last : = new ; 

end; 

end 

else begin 

subregion( . last. ). active := true; 
joined := false ; 
last := new ; 

end; 

end 

else begin 

subregion(. last. ). active := true; 
joined : = false ; 
last : = new ; 

end; 

end; 

end; (* column combine finished *) 
tabc_dif ference := tabc_dif f erence + increment; 
end; (* while *) 

end; (* procedure combine *) 

procedure cal_avg_dev( var subregion: subregion_array; 

var points: point_array; var avg:real); 

var 

sum: real; 
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i, j ,gp: integer; 
begin 

sum : = 0. 0; 

for i: = 2 to row-1 do 

for j:= 2 to col-1 do 

begin 

gp : = points(. i, j. ). gp; 

sum := sura + abs(( subregion( . gp. ) • 

subregion(. gp. ). b*j + 
subregion(. gp. ). c) 

- points(. i, j. ). el); 



end; 

avg := sum / ((row-2)*( col-2)); 
writeln( * average deviation = *,avg : 15); 
end; (* sub procedure cal_average_deviation *) 



procedure writedata(points: point_array); 
var 

r,c : integer; 
tempa , tempb, tempc : real; 
begin 

for r: = 2 to row-1 do 
for c: = 2 to col-1 do 
begin 

tempa := subregion( . points( . r ,c. ). gp. ). a; 

tempb := subregion( . points( . r , c. ). gp. ) . b; 

tempc := subregion( . points( . r , c. ) . gp. ) . c; 

points( . r , c. ). el := trunc (tempa*r + tempb*c + tempc) 

end; 

for c : = 1 to col do 
for r : = row downto 1 do 

writeln(doutl,points( . r,c. ). el); 

end; 

(* procedure writedata ’^) 



procedure writetempdata(points: point_array) ; 
var 

k,r,c,rl,r2,cl,c2 : integer; 

tempa, tempb, tempc : real; 
current : nodepointer; 
begin 

for r: = 2 to row-1 do 
for c: = 2 to col-1 do 

points(. r,c. ). el := 1 ; 
for k: = 1 to counter do 
if subregion( . k. ). active then 

begin 

current := subregion( . k. ). list ; 
while current <> nil do begin 
rl: =current“*. r; 
cl: =current-*. c; 

if current"’, next <> nil then begin 
r2: =current"*. next-*, r; 
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c2: =current-', next^. c; 

end 

else begin 
r2 := rl; 
c2 : = cl; 
end; 

if rl > r2 then begin r: = r2; r2: = rl; rl:= r; 



if cl > c2 then begin c: = c2; c2: = cl; cl:= c; 

tempa := subregion( . k. ) . a; 
tempb := subregion( . k. ) . b; 
tempc := subregion( . k. ) . c; 
if rl=r2 then for c : = cl to c2 do 
pointsC . rl,c. )• •= trunc (tempa*rl + tempb*c 
if cl=c2 then for r := rl to r2 do 
pointsC . r ,cl. )• trunc (tempa*r + tempb*cl 
current ;= current-’, next; 
end; 
end; 

for c : = 1 to col do 
for r : = row downto 1 do 

writeln(dout2,points( . r,c. ). el); 

end; 

(* procedure writetempdata *) 



procedure check; 
var 

num,r,c,i : integer; 
begin 

for c: = 1 to counter do temparr(.c. ) := false; (* 
num : = 0; 

for i : = 1 to counter do 

if subregion(. i. ). active then 
begin 

r:= i ; 

temparr(.r. ) := true ; 
num : = num -l-l ; 
while r <> 0 do 
begin 

with subregion(. r. ) do 
temparr(.r. ) := true; 
r: = subregion( . r. ). side ; 

end; 

end; 

for c: = 1 to counter do temparr(.c. ) := false; (* 
num := 0; 

for i : = 1 to counter do 

if subregion(. i. ). active then 
begin 

r:= i ; 

temparr(.r. ) := true ; 
num : = num +1 ; 



end; 

end; 

+ tempc); 
+ tempc); 



check *) 



check *) 
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while r <> 0 do 
begin 

with subregion( . r. ) do 
temparr(.r. ) := true; 
r: = subregion( . r. )• side ; 

end; 

end; 

writeln(*num of gps = counter :3,* lowbound=* , low_bound: 3 , 

* upper_bound=* ,upper_bound: 3, * * max_sd=* ,max_sd: 3 ) 

writeln( * unexamined gps are =>*); 
c: = 0; 

for i : = 1 to counter do if temparr(. i. ) = true then c: = c+1 
else write( ,i: 3, *** ); 
writeln; 

writeln(*# of gps examined =>’,c:4,* * abc_difference=* , 
abc_dif f erence) ; 

writeln(*# of combined gps =>’ ,num :4); 
end; (* procedure check *) 
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( ^= ==^-= =====:r = ==;=rr=:^rT==r=== = =-^— r=====^ = 

Main Program 

begin 

page; 

initialize ; 
new( top); 
top-*, up : = nil ; 

reset(data, Mata ell*); 
reset(datam, *data magi*); 

rewrite(doutl , * data outl*); data after combine process *) 

rewrite( dout2 , * tdata outl*); (* data of boundary *) 

for c ; = 1 to col do 
for r : = row downto 1 do 

readln(data,points( . r. )(. c. ). el); 
for c : = 1 to col do 
for r : = row downto 1 do 

readlnC dat am, point s( . r. )( • c. )• ni^g)^ 

makesub(2,(row-l) ,2,(col-l) ,top); 

printsubl(points); (* print the subregions created by makesub *) 

for i : = 1 to counter do 
with subregion(. i. ) do 
begin 

get const (points ,cr 1 , crh, ccl , ccr ,ok, ee, ex,ey); 
gauss(ok,a,b,c); 

sd := sqrt(variance(a,b,c,ee,ex,ey,ok)); 

(* sd for region i after division *) 



constarr : = ok ; 
end; 

comb ine( points , subregion); 
get_nodes( subregion) ; 

printsub2(points); (* print the subregions created by combine ’^^) 
wr it edata( points ) ; 

(* write modified elevation data into file *data out//’*) 
wr it etempdata( points); 

(* write elevation data into file ’tempdata out//’ *) 
cal_avg_dev( subregion , points , avg) ; 

(* check; *) 

end. 
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APPENDIX B. PASCAL CODE OF STRICT BOTTOM-UP ALGORITHM 



(*$s450000*) 



Author 

Date 

Input 



Output = 



Computer 

Descript 



= Yee, Seung Hee. 

= 3 June 1988. 

= 1. Elevation data acquired from Fort. Hunter Liggett 
data base (40 by 40 square). 

2. file of magnitude of gradient of the above elevation 
data created by using the program in Appendix 3. 

3. file of curvature of points created by using the 
program in Appendix 3. 

1. Elevation data of the model created by this program. 

2. Elevation data of the boundary points, other points 
are given the value one so that when we display 
this model in three-dimensional pictures, we can 
only see the boundaries of the patches. 

= Waterloo Pascal on IBM 370/3033AP, 

Naval Postgraduate School, 
ion = 

This program is the Pascal implementation of the 
Joint top-down and bottom-up terrain modeling. 






PROGRAM Bottom_up( input , output ) ; 
const 

low_bound = 0. 1; magnitude < 10 % then level slope *) 

(* magnitude > 50 % then unsafe slope *) 
up_bound = 0. 6; (* else safe slope *) 

magnitude_ratio = 0. 2; (* ratio of magnitude : 

dif of two mags / old mag *) 

abc_comparison_const =20; 
max_group_num = 400; 
row = 40; 
col = 40; 
type 

elrange = 1. . 10000; 

mag^type = ( level, safe, unsafe); 

point = record 

r ,c: integer; 
end; 

nodepointer = -•node; 
node = record 

r,c : integer; 
next : nodepointer; 
end; 

point_rec = record 

el : elrange ; 
mag : real; 
magtype: mag_type; 
cur: char; 

gp : 0. . max_group_num; 
ptnext : point; 
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end; 

point^array = array ( . 1, , row, 1. . col. ) of point_rec; 
polynomial_array = array ( . 1. . 3 , 1. . 4. ) of real; 
coef_record = record 

start : point; 

upper_left : point; 

a,b,c : real ; 

sd : real ; 

ee,ex,ey : real; 

side : 0. . max_group_num; 

poly_array : polynomial_array; 

active : boolean; 

list : nodepointer ; 

numpts : integer; 

avgraag: real; 

end; 

coef_array = array ( . 1. . max_group_num. ) of coef_record ; 

var 

points : point_array ; 
subregion : coef_array ; 
ok : polynoraial_array; 

data,datam,datac,doutl ,dout2 ,dout3 ,dout4 : text; 
counter , r, c, i ,num ; integer; 
tempmag : real; 

function variance (a,b,c,ee,ex,ey: real; ok: polynomial_array): real; 
begin 

variance := ( ee-( 2*a*ex) -( 2’^b*ey) -( 2*c*ok( . 3 ,4. )) 

+( a*a*ok( . 1 , 1. ) )+( 2*a*b*ok( .2,1.)) 
+(2’Vb*c*ok(. 2,3. ))+(2*a*c*ok(. 1,3. )) 
+(b*b*ok(. 2,2. ))+(c*c*ok(. 3,3. ))) 
/(ok(.3,3. )-l); 

end; 

procedure initialize; 
var 

r,c .-integer; 
begin 

counter : = 0 ; 

for c : = 1 to col do 

for r : = 1 to row do 

begin 

points(. r,c. ). gp : = 0 ; 
points(, r,c. ). ptnext. r : = 0 ; 
points( . r,c. ). ptnext. c : = 0 ; 
end; 

for r : = 1 to max_group_num do 
begin 

subregion( . r. ). a := -999; 
subregion( . r. ). b := -999; 
subregion( . r. ). c := -999; 
subregionC . r. ). active := true; 
subregion( . r. ). side := 0; 
subregionC . r. ). sd := -999; 

end; 

end; (* procedure initialize *) 
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procedure 

begin 



end; 



getconst( r,c: integer; 

var arr: polynomial^array; var ee,ex,ey : real); 



arr( 


1,1. 


) 


=arr( , 


, 1,1. 


)+(r*r); 




arr( 


1,2. 


) 


=arr( . 


, 1,2. 


)+(c*r); 




arr( 


1,3. 


) 


=arr( . 


1,3. 


)+ r; 




arr( 


1,4. 


) 


=arr( , 


1,4. 


)+(points( . r ,c. 


).el*r); 


arr( 


2,1. 


) 


=arr( . 


2,1. 


)+(r*c); 




arr( . 


2,2. 


) 


=arr( . 


2,2. 


)+(c*c); 




arr( . 


2,3. 


) 


=arr( . 


2,3. 


)+ c; 




arr( . 


2,4. 


) 


=arr( . 


2,4. 


)+(points(. r,c. 


). el*c); 


arr( . 


3,1. 


) 


=arr( . 


3,1. 


)+r; 




arr( . 


3,2. 


) 


=arr( , 


3,2. 


)+c; 




arr( . 


3,3. 


) 


=arr( . 


3,3. 


)+i; 




arr( . 


3,4. 


) 


=arr( , 


3,4. 


)+points(. r ,c. ) 


.el; 


ee : = 


= ee 


+ 


sqr(points(. r,c. ). el); 




ex : = 


= ex 


+ 


(points( . r 


,c. ). el * r ); 




ey : = 


= ey 


+ 


(points( . r 


,c. ). el * c ); 





(* procedure get const *) 



procedure gauss(ok: polynomial_array; var aa,bb,cc: real); 



var 

ij : integer; 

temp : array (.1..4. ) of real; 

t : real; error checking purpose *) 

begin 

i := 1; 

while ok(.i,l.) = 0. 0 do i := i + 1 ; (* row rotation *) 

for j := 1 to 4 do 

begin 

temp(. j. ) : = ok(. 1, j. ); 
ok(. 1, j. ) : = ok(. i, j. ); 
ok(. i^ j. ) : = temp(. j. ); 
end; 

t : = ok(. 1,4. ); 

for i : = 2 to 3 do 

for j : = 4 downto 1 do 

ok(.i,j.) := ok(. i, j. )*ok(. 1, 1. ) + ok( . 1 , j. )*( -ok( . i, 1. ) ); 
i : = 2; 

while ok(.i,2. ) = 0.0 do i:=i+l; 

for j := 1 to 4 do 

begin 

temp(. j. ) := ok(. 2, j. ); 
ok(. 2, j. ) := ok(. i,j. ); 
ok(. i, j. ) : = temp(. j. ); 
end; 

for j : = 4 downto 1 do 
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ok(.3,j.) := ok(.3,j.)*ok(.2,2. ) + ok( . 2 , j . -ok( . 3 , 2. ) ) ; 



if ok(.3,3. ) <> 0.0 then 
begin 

ok(.334.) := ok(.3,4.) / ok(.3,3.); 

ok(.2,4. ) := ( ok(.2,4. ) - ok( . 3 ,4. )*ok( . 2, 3. ) ) / ok(.2,2. ) 
ok(.l,4. ) := ( ok(.l,4. ) - ok( . 2 ,4. )*ok( . 1 , 2. ) 

- ok(. 3,4. )*ok(. 1,3. )) / ok(.l,l. ) ; 
aa := ok(.l,4. ); (* solution *) 
bb := ok(.2,4. ); (* solution *) 
cc := ok(.3,4. ); solution *) 

(*writeln( ’ gauss error , 

t-( aa*ok(. 1,1. )+bb*ok(. 1,2. )+ cc*ok(. 1,3. )))*) 

end; 

end; (* procedure gauss *) 

procedure prefix(nura: integer); 
begin 

case ((num div 100) mod 10) of 

0 : write( * * ); 

1 : write( \ ^ 

2 : write( * : ^ ); 

3 : write( ^ - ^ ); 

4 : write(*+^); 

5 : write(’%’); 

6 : write(^**); 

7 : write( ); 

8 : write( ^ ? * ); 

9 : write( ^ = * ); 

end; 

end; (* procedure prefix *) 

procedure printsub(points: point_array); 
var 

r,c : integer; 
begin 
page; 

for r : = 1 to row do 
for c : = 1 to col do 
begin 

if c<> col then begin 

pref ix( points ( . r,c. ). gp); 
write( (points( . r , c. ) . gp mod 100): 2); 
end 

else begin 

pref ix( points ( . r,c. ). gp); 
writeln((points(. r,c. ). gp mod 100): 2); 
writeln; 
end; 

end; 

end; procedure printsub *) 



procedure pr intsub2( points : point_array ) ; 
var 

k,r ,c,rl,r2,cl,c2 : integer; 
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current : nodepointer; 
begin 
page; 

for r: = 2 to row-1 do 
for c:= 2 to col-1 do 

points (. r,c. )• gp : = 0 ; 
for k: = 1 to counter do 
if subregion( . k. )• active then 
begin 

current := subregion( . k. )• list ; 

while current <> nil do 

begin 

rl: =current“*, r; 
cl: =current-’. c; 

if current-*, next <> nil then begin 
r2: =current“*. next-*, r; 
c2: =current-*. next-*, c; 

end 

else begin 
r2 := rl; 
c2 := cl; 
end; 

if rl > r2 then begin r: = r2; r2:= rl; rl:= r; 
if cl > c2 then begin c: = c2; c2: = cl; cl: = c; 

if rl=r2 then for c : = cl to c2 do 
points( . rl j c. ). gp : = k; 
if cl=c2 then for r : = rl to r2 do 
points( . r , cl. ). gp : = k; 
current := current-*, next; 
end; 
end; 

for r : = 1 to row do 
for c : = 1 to col do 
begin 

if c<> col then begin 

prefix(points(. r,c. ). gp); 
write((points(. r ,c. ). gp mod 100): 2); 
end 

else begin 

prefix( points ( . r , c. ). gp); 
writeln((points(. r,c. ). gp mod 100): 2); 
writeln; 
end; 

end; 

end; 

(* procedure printsub2 *) 

procedure printmagcur; 
var 

r,c : integer; 

begin 
(* page*) 

for r : = 1 to row do 
for c : = 1 to col do 
begin 



end; 

end; 
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if c<> col then begin 

case points(. r ,c. )• nisgtype of 

level : write( ’ 1* ,points(. r,c. ). cur, * * ); 
safe : write( * s * ,i^oints( . r , c. ). cur , * * ) ; 
unsafe : write( ’u’ ,points(. r ,c. ). cur, * * ); 
end; (* case *) 
end 

else begin 

case points( . r ,c. )• inagtype of 

level : write ln( ' 1* ,points(. r,c. ). cur, * * ); 

safe : writeln( * s* ,points(. r ,c. ). cur, * ’ ); 

unsafe; write ln( *u* ,points(. r,c. ). cur, * * ); 

end; 

writeln; 

end; 

end; 

end; (* procedure printmagcur *) 

procedure makeone_rowl( var pointl ,point2 : point); 

begin 

points( . pointl. r, pointl. c. ). ptnext. r : = point2. r; 
points( . point 1. r , point 1. c. ). ptnext. c := point2. c; 
points (. point 2. r, point 2. c. ). gp : = 

point s( . pointl. r, pointl. c. ). gp; 
with subregion( . counter. ) do 
numpts ; = numpts +1; 



end; (* sub procedure make one_l *) 



procedure makeone_row2( var point 1 ,point2 ; point); 
begin 

points( . pointl. r , point 1. c. ). ptnext. r : = point2. r; 
points( . pointl. r, pointl. c. ). ptnext. c := point2. c; 
points( . point2. r ,point2. c. ). gp : = 

points( . pointl. r, pointl. c. ). gp; 
with subregion( . counter. ) do 
begin 

numpts : = numpts +1; 

avgmag := avgmag*(numpts-l)/numpts + 

points( . point2. r ,point2. c. ). mag/numpts; 

end; 

end; (* sub procedure make one_2 *) 



procedure rowlink; 
var 

i,j : integer; 

point l,point2 : point; 
begin (* sub procedure rowlink *) 

counter ; = 0; 
for i : = 2 to (row-1) do 
begin 
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counter : = counter + 1 ; 
points( . i , 2. )• gP counter; 

subregion(. counter. ). start. r : = i; 
subregion( . counter. ). start, c := 2; 
subregion( . counter. ) . upper_lef t : = 

subregion(. counter. ). start ; 
subregion( . counter. ). avgraag : = points( . i,2. ). mag; 
subregion( . counter. ). numpts := 1 ; 
for j : = 2 to (col-2) do 
begin 

pointl. r := i; 
pointl. c : = j ; 
point2. r : = i; 
point2. c : = j+1; 

if ((points(. i, j+1. ). magtype <> safe) and 

( points (. i,j. ). magtype = points(. i, j + 1. ). magtype)) then 
makeone_rowl( point 1 , point 2) 

(* (level level) or (unsafe unsafe) *) 
else if ( (points( . i,j. ). magtype = safe) and 

(points(. i, j+1. ). magtype = safe)) and 
(points(. i, j. ). cur = points(. i, j+1. ). cur) then 
if (abs(points(. i, j+1. ). mag- 

subregion(. counter. ). avgmag)/ 
subregion(. counter. ). avgmag) < 
magnitude_ratio then 
makeone_row2 ( point 1 , point 2 ) 
else begin 

counter : = counter +1 ; 
points(. i, j+1. ). gp := counter; 
subregion( . counter. ). start := point2; 
subregion( . counter. ). upper_left := point2 ; 
subregion(. counter. ). avgmag : = 

points(. i, j+1. ). mag; 
subregion( . counter. ) . numpts : = 1 ; 
pointl := point2; 

end 

else begin 

counter : = counter +1 ; 
points( . i, j+1. ). gp := counter; 
subregion( . counter. ). start := point2; 
subregion( . counter. ). upper_left := point2 ; 
subregion(. counter. ). avgmag : = 

points(. i,j+l. ).mag; 
subregion( . counter. ). numpts := 1 ; 
pointl := point2; 

end; 

end; 

end; 

end; (* sub procedure row link *) 



procedure makeone_col( last , new : point); 
var 

ii,jj : integer; 
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tail, current : point; 



begin 

subregion(. points(. new. r, new. c. )• gP- )• active := false ; 
with subregion(. points( . last. r, last. c. ). gp. ) do 
begin 

numpts: =numpts + 

subregion(. points(. new. r,new. c. ). gp. ). numpts; 
avgmag: =avgmag*( numpts - 

subregion(. points (. new. r,new. c. ). gp. ). numpts)/ 
numpts + 

subregion( . points( . new. r ,new. c. ). gp. ). avgmag* 
subregion(. points(. new. r,new. c. ). gp. ). numpts /numpts; 
end; (* reassigning number of points and average magnitude *) 

if (subregion(. points(. new. r,new. c. ). gp. ). upper_left. r + 
subregion( . points( . new. r ,new. c. ). gp. ). upper_left. c) < 

(subregion( . points(. last, r, last. c. ). gp. ). upper_left. r + 
subregion(. points(. last. r,last. c. ). gp. ). upper_left. c) then 
subregion(. points(. last. r, last. c. ). gp. ). upper_left : = 
subregion(. points(. new. r,new. c. ). gp. ). upper_left ; 

current := points( . last. r , last. c. ). ptnext ; 

if current. r = 0 then (* or current. c = 0 *) 

begin 

points( . last, r , last. c. ) . ptnext : = 

subregion(. points(. new. r,new. c. ). gp. ). start ; 
current : = subregion(. points(. new. r,new. c. ). gp. ). start ; 

end 

else begin 

while current. r <> 0 do 
begin 

tail := current; 

current := points(. current, r, current, c. ). ptnext; 
end; (* traverse ptnext of lastpt link 

to the end of link *) 

points( . tail, r , tail. c. ). ptnext: = (* now link the tail 

points to the head of*) 
subregion(. points(. new. r,new. c. ). gp. ). start ; 

(* newpt link *) 

current := (* into existing link *) 

subregion(. points(. new. r,new. c. ). gp. ). start ; 

end; 

while current, r <> 0 do (* reassigning gp # to new pts *) 
begin 

points(. current, r, current, c. ). gp : = points(. last. r , last. c. ). gp; 
current : = points ( . current, r , current, c. ) . ptnext; 
end; (* while *) 

end; (* procedure makeone_col *) 



procedure col link; 
var 

i,j : integer; 
point 1 ,point2 : point; 
begin (* sub procedure col link *) 
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for j : = 2 to (col-1) do 
begin 

for i : = 2 to (row-2) do 
begin 

point 1. r : = i; 
pointl. c : = j; 
point2. r := i+1; 
point2. c : = j; 

if (points(. i+1, j. )• gP points( . i , j . ) • gp) then 

if ( (points (. i+1 , j .)• niagtype <> safe) and 

( point s (. i,j. )• niagtype = points( . i+1 , j . ) • roagtype) ) then 
makeone_col( pointl , point 2) 

(* ( level, level) or (unsafe, unsafe) *) 
else if ((points(. i, j. )• loagtype = safe) and 

(points(. i+1, j. )• magtype = safe)) then 
if (points(. i, j. )• cur = points( . i+1 , j . ) . cur) then 
if abs(subregion( . points(. i+1, j, ), gp. ). avgmag- 
subregion(. points(. i, j. )• gp* )• avgmag)/ 
subr egion( . points ( . i , j . ) . gp. ) . avgmag < 
magnitude_r at io 

then 

makeone_col( point 1 , point 2 ) ; 
pointl := point2; 
end; 
end; 

end; (* sub procedure collink *) 
procedure do_combine( last , new : point); 
var 

ii,jj : integer; 
tail, current : point; 

begin 

subregion( . points( . new. r , new, c. ). gp. ). active := false ; 
with subregion( . points( . last. r, last. c. ). gp. ) do 
begin 

nump t s : =nump t s + 

subregion(. points(. new. r,new. c. ). gp. ). numpts; 
end; (* reassigning number of points *) 

if (subregion( . points( . new. r ,new. c. ). gp. ). upper_left. r + 
subregion(. points(. new. r ,new. c. ), gp. ). upper_left. c) < 
(subregion(. points(. last. r,last. c. ). gp. ). upper_left. r + 
subregion(. points (. last. r,last. c. ). gp. ). upper_left. c) then 
subregion( . points( . last. r , last. c. ) . gp. ) . upper_left : = 
subregion(. points(. new. r,new. c. ). gp. ). upper_left ; 



current := points(. last. r, last. c. ). ptnext ; 

if current. r = 0 then (* or current. c = 0 *) 

begin 

points( . last. r , last. c. ) . ptnext : = 

subregion( . points( . new. r ,new. c. ). gp, ). start ; 
current : = subregion(. points(. new. r,new. c. ). gp. ). start ; 

end 
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else begin 

while current, r <> 0 do 
begin 

tail := current; 

current := points( . current, r, current, c. ). ptnext; 
end; (* traverse ptnext to link the head *) 

points( . tail, retail, c. ). ptnext: = (* of the *) 

subregion( . points( . new. r ,new. c. ). gp. ). start ; 

(* new point *) 

current := (* into existing link *) 

subregion( . points ( . new. r ,new. c. ). gp. ). start ; 

end; 

while current, r <> 0 do (* reassigning gp # to new pts *) 
begin 

points(. current, r, current, c. ). gp : = points(. last. r,last. c. ). gp; 
current := points( . current, r , current, c. ). ptnext; 
end; (* while 



for ii : = 1 to 3 do 
for j j : = 1 to 4 do 

subregion( . points( . last. r , last. c. ) . gp 
subregion( . points( . last. r , last. c 
subregion( . points ( . new. r ,new. c. ) 
subregionC . points( . last, r , last. c. ). gp 
subregion( . points( . last. r , last. c 
subregionC . points( . new. r ,new. c. ) 
subregionC . pointsC . last, r , last. c. ). gp 
subregionC • pointsC . last. r , last. c 
subregionC . pointsC . new. r ,new. c. ) 
subregionC • pointsC • last, r , last. c. ). gp 
subregionC . pointsC . last. r , last. c 
subregionC • pointsC • new. r ,new. c. ) 
with subregionC • pointsC • last. r , last. c 
begin 



). poly_array(. ii, j j. ) : = 

). gp. )• poly_array(. ii, j j. )+ 
gp. ). poly_array(. ii, j j. ); 

) . ee : = 

). gp. ). ee + 
gp. ). ee ; 

) . ex : = 

). gp. ). ex + 
gp. ). ex ; 

) . ey : = 

). gp. ). ey + 
gp- )• ey ; 

). gp. ) do 



sd := sqrtCvarianceCa,b,c,ee,ex,ey,poly_array)); 

end; 

end; C’"^ procedure do_combine *) 



procedure combineCvar pt: point_array; 

var subregion: coef_array); 

var 

i,j : integer; 
last, new : point; 
joinabCjjoinsd, joined : boolean; 

procedure compareabcC var join: boolean); 
var 

temp ; real ; 
begin 

temp ;= sqrtC sqrC subregionC . pointsC . new. r , new. c. ). gp. ). a - 
subregionC. pointsC. last. r,last. c. ). gp. ). a) + 
sqrC subregionC . pointsC . new. r ,new. c. ). gp. ). b - 
subregionC • pointsC . last, r , last. c. ). gp. ). b) 4* 
sqrC subregionC • pointsC • new. r ,new. c. ). gp. ). c - 
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subregion( . point s( . last. r,last. c. ). gp. ) . c) ); 
if temp <= abc_comparison_const then join := true 
else join := false; 

end; sub procedure compareabc 



begin (* procedure combine *) 

joined : = false; 
for i : = 2 to (row-2) do 
begin 

last, r : = i; 
last, c : = 2; 

for j : = 2 to (col-2) do 
begin 

new. r : = i; 
new. c : = j + 1; 

if ((i=row-2) and (j=col-2)) and (not joined) then 

subregion(. points (. new. r ,new. c. ). gp. ). active : = 

true; 

if (points( . new. r ,new. c. ). gp <> points( . last. r , last. c. ). gp) 
then begin 

compareabc( joinabc); 
if joinabc 
then begin 

do_combine( last ,new); 
joined := true; 
end 

else begin 

subregion(. points(. last, r , last. c. ). gp. ). active : = 

true; 

joined : = false ; 
last : = new ; 
end; (* joinabc *) 

end 

else begin 

subregion( . points( . last. r , last. c. ). gp. ). active : = 

true; 



joined := false ; 
last : = new ; 

end; (* last <> new *) 
end; (* end j *) 

end; (* end i *)(* row combine finished *) 

for j : = 2 to (row-2) do 

begin 

last, r : = 2; 
last. c : = j; 

for i : = 2 to (col-2) do 
begin 

new. r : = i+1; 
new. c : = j ; 

if (points( . new. r ,new. c. ). gp <> points( . last, r, last. c. ). gp) 
then begin 

compareabc( j oinabc) ; 
if joinabc 
then begin 
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do_combine( last , new); 
joined := true; 

end 

else begin 

subregion( . points( . last. r , last. c. ). gp. ). active 

true; 

joined : = false ; 
last := new ; 
end; (* join abc *) 

end 

else begin 

subregion(. points(. last. r,last. c. ). gp. ). active : = 

true; 

joined : = false ; 
last : = new ; 

end; (* last <> new *) 

end; 

end; (* column combine finished *) 
end; (* procedure combine *) 



procedure get_boundar ies ( var subregion; coef_ar ray ) ; 
type 

direction = ( up, down, right , left ); 

var 

i,j,k,r,c : integer; 

current , head, p : nodepointer; 
temp : direction ; 



procedure get_next(var p : nodepointer; 

var temp : direction); 



var 



i,j ; integer; 
done : boolean ; 



procedure do_dir(var done: boolean 
begin 

new(p); 

P-. r : = i ; 
p-'-c := j ; 
done : = true; 
case dir of 

temp : = up; 
temp ; = down; 
temp := right; 
temp : = left ; 

(* case *) 

sub_sub_ procedure do_dir 



dir: direction ); 



up 

down 
right 
left 
end; 
end; (* 



begin (* sub procedure get next *) 
done := false ; 
i := P-. r; 
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j := p-’- c; 
case temp 
right : 



left : 



up 



down : 



of 

begin 

j •= j + 1; 

while not done do begin 

if points(. i-1, j. )• gp = k 

then do_dir(done,up) (* go up and done *) 
else if( ( subregion( . k. ). list-*, r = i) 

and (subregion( . k. )• list-*, c = j)) then 
begin (* return to org point *) 

new(p); 

P-. r := i ; 

:= j ; 

done : = true; 



end 

else if (points(. i, j+1. )• gp = k) 

then j:= j+1 (* go right *) 

else if (points( . i+1, j. )• gp=k ) 

then do_dir( done, down) (* go down and done *) 
else do_dir(done, left); (* go back and done *) 

end; (* while *) 
end; 
begin 

j •= j-i; 

while not done do begin 

if points( . i+1 , j. ). gp = k 

then do_dir( done, down) (* go down and done*) 
else if (points(. i, j-1. )• gp = k) 



then j : = j-1 
else if (points(. i-1, j. )• gp = k) 



(* go left *) 



(* go up and done *) 
(* go back and done *) 



then do_dir(done,up) 
else do_dir( done, right); 
end; (* while *) 
end; 
begin 

i := i-1; 

while not done do begin 

if points(. i, j-1. )• gp = k 

then do_dir( done, left) (* go left and done *) 
else if (points(. i-1, j. ). gp = k) 

then i := i-1 (* go up *) 

else if (pointsC . i, j+1. ). gp = k ) 

then do_dir(done, right )(* go right and done*) 
else do_dir( done, down); (* go down and done *) 

end; (* while *) 
end; 
begin 

i := i+1; 

while not done do begin 

if pointsC . i, j+1. ). gp = k 

then do_dir (done, right )(* go right and done*) 
else if pointsC . i+l,j. ). gp = 



then i : = i+1 



(* go down *) 
= k 



else if pointsC . i,j -1. ). gp = 

then do_dir(done, left )(* go left and done*) 
else do_dir(done,up); (* go up and done *) 
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end; (* while *) 
end; 

end; (* case *) 

end; (* sub procedure get_next *) 



begin procedure get boundaries *) 

for k : = 1 to counter do 
if subregion( . k. )• active then 
begin 

new(p); 

p-'. r : =subregion( . k. )• upper^left. r ; 
p-'. c : =subregion(. k. )• upper_left. c ; 
while ( pointsC . p-'. r-1 ,p-’. c. )• gp = k) do 

p-’. r := p*’. r-1; (* move to the upper most row of the gp *) 

while ( points( . p*'. r ,p*'. c-1. ) • gp == k) do 

p-’. c := p*’. c-1; (* move to the left most col of the gp *) 

subregion( . kO- list : = p ; 

if points(. p-*. r ,p-*. c+1. )• gp = k then 

begin (* starting toward right groups *) 

current : = p ; 
temp : = right ; 
repeat 

ge t_next ( p , t emp ) ; 
current-*, next : = p ; 
current : = p ; 

until (subregion(. k. ). list-’, r = p“'. r) 

and (subregion(. k. ). list-’, c = p-’.c); 
if (temp=down) and (points( . p-’. r+1 ,p-*. c. ). gp=k) then 
begin (* right and down points gp *) 

current := p ; (* hinged on start point *) 

temp : = down ; 
repeat 

get_next ( p , temp) ; 
current-’, next : = p ; 
current : = p ; 

until (subregion(. k. ). list-’, r = p"*. r) 

and (subregionC . k. ). list-*, c = p-*. c); 
p-*. next : = nil; 

end 

else 

p-*. next := nil; (* right but not down gp *) 



end 

else if points(. p-*. r+ljp-*. c. ). gp = k then 
begin starting toward down group *) 

current : = p ; 
temp : = down ; 
repeat 

get^next ( p , temp ) ; 
current-*, next : = p ; 
current : = p ; 

until ( subregionC . k. ). list-*, r = p"'. r) 

and (subregion(. k. ). list-*, c = p*’. c); 
p-’. next : = nil; 
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end 

else 

p-*.next := nil; (* one point gp *) 

end; 



the following segment of code may be used for checking purpose 
for k : = 1 to counter do 
if subregion( . k. )• active then 
begin 

write( ' Boundary of group *,k:3/is =>*); 
current := subregion( . k. ) • list; 
while current <> nil do 
begin 

wr it e( current*^, r: 2, * , * , current-*, c: 2/=>* ); 
current := current-*, next; 
end; 

writeln; 

end; 



end; (* procedure get_boundaries *) 



.V) 



procedure calc_coef fs; 
var 

k, i, j ,r jCjgpnumjptnum : integer; 
t,tl : point; 

rowchanged, colchanged: boolean; 
begin 

gpnum : = 0; 
ptnura : = 0; 

for k : = 1 to counter do 

if subregion(. k. ). active then 
begin 

rowchanged := false; 
colchanged := false; 
for i : = 1 to 3 do 
for j := 1 to 4 do 

subregion(. k. ). poly_array(. i, j. ) := 0; 
subregionC . k. ). ee := 0; 
subregion( . k. ). ex := 0; 
subregionC . k. ). ey := 0; 
ptnum := ptnura +subregion( . k. ) . nurapts ; 

(* num of linked pt *) 

gpnura := gpnum +1 ; (* num of active gp *) 

t := subregionC . k. ). start; 
while t. r <> 0 do 
begin 

tl := t; 

with subregionC . k. ) do 

getconstC t. r, t. c,poly_array,ee,ex,ey); 
t:= pointsC . t. r j t. c. ). ptnext ; 
if Ct.roO) and Ct.r<>tl. r) then 
rowchanged : = true; 
if Ct.coO) and Ct.c<>tl. c) then 
colchanged : = true; 

end; 

with subregionC . k. ) do 
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begin 

if rowchanged and colchanged then 

gauss(poly_array ,a,b,c) (* normal points *) 

else if not( rowchanged) and not ( colchanged) then 
begin 

a : = 0; 
b : = 0; 

c := poly_array(. 3,4. ); 

end (* one point gp *) 

else if not( rowchanged) and colchanged then 
begin (* horizontal line group *) 

a ; = 0; 

for j:= 4 downto 2 do 
poly_array(. 3, j. ) : = 

poly_array(. 3, j. ) * poly_array( . 2 , 2. ) + 
poly_array(. 2, j. ) * (-poly_array(. 3,2. )); 
poly_array(. 3,4. ) : = 

poly_array(. 3,4. ) / 
poly_array(. 3,3. ); 

poly_array(. 2,4. ) := (poly_array(. 2,4. ) 

poly_array(. 3,4. ) * 
poly^arrayO 2,3. )) / 

poly_array(. 2,2. ); 
b := poly_array(. 2,4. ); 
c ;= poly_array(. 3,4. ); 

end 

else if rowchanged and not( colchanged) then 
begin (* vertical line group *) 

b : = 0; 

for j:= 4 downto 1 do 
poly_array(. 3, j. ) : = 

poly_array(. 3, j. )* 
poly_array( .1,1.) + 
poly_array(. 1, j. )* 

(-poly_array(. 3,1. )); 
poly_array(. 3,4. ) : = 

poly_array(. 3,4. ) / 
poly_array(. 3,3. ); 

poly_array( . 1 ,4. ) := (poly_array(. 1 ,4. ) - 
poly_array(. 3,4. )* 
poly^arrayO 1,3. ))/ 
poly_array(. 1,1.); 
a := poly_array( . 1 ,4. ); 
c := poly_array( . 3 ,4. ); 
end; 

(* writeln( * a , b , c , = * , a: 10 , ' *,b:10,* ’ ,c: 10) *) 

end; (* with subregion( . k. ) *) 

end; (* if subregion( . k. ). active *) 



writelnC * active pts=* ,ptnura: 5 , * * active gpnum =*,gpnum:3); 
end; (* procedure calc_coeffs *) 



procedure write_planar_data(points: point_array); 
var 

r,c ; integer; 
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tempa,tempb, tempo : real; 
begin 

for r: = 2 to row-1 do 
for c: = 2 to col-1 do 
begin 

tempa := subregion( . points(. r ,c. )• gP« 

tempb := subregion( . points( . r , c. ) • gP- ) • b; 
tempo := subregion( . points( . r , c. ) • gP- ) • c; 
points( . r ,0. )• trunc (tempa*r + tempb*c + tempo); 

end; 

for c : = 1 to col do 
for r : = row downto 1 do 

writeln(doutl ,points( . r,c. ). el); 

end; 

(* procedure write_planar_data *) 



procedure write_boundary_data( points: point_array); 
var 

k,r,c,rl,r 2 ,cl 5 c 2 : integer; 

tempa, tempb, tempo : real; 
current : nodepointer; 
begin 

for r:= 2 to row-1 do 
for c: = 2 to col-1 do 

points(. r,c. )• el := 1 ; 
for k: = 1 to counter do 
if subregion( . k. )• active then 
begin 

current := subregion( . k. )• list ; 

while current <> nil do 

begin 

rl: =current“'. r; 
cl: =current"'. c; 

if current-*, next <> nil then begin 
r2: =current-*. next-*, r; 
c2: =current-'. next-*, c; 

end 

else begin 
r2 := rl; 
c2 : = cl; 
end; 

if rl > r2 then begin r: = r2; r2:= rl; rl:= r; end; 
if cl > c2 then begin c: = c2; c2:= cl; cl:= c; end; 

tempa := subregion( . k. ). a; 
tempb := subregion( . k. ) . b; 
tempo := subregion( . k. ). c; 
if rl=r2 then for c : = cl to c2 do 

points( . rl , c. ). el := trunc (tempa*rl + tempb*c + tempo) 
if cl=c2 then for r := rl to r2 do 

points( . r ,cl. ). el := trunc (tempa*r + tempb*cl + tempo) 
current := current-*, next; 
end; 
end; 

for c : = 1 to col do 



82 



for r : = row downto 1 do 

writeln(dout2 , points ( . r , c. ) . el); 
end; (* write_boundary_data *) 

(* procedure write_boundary_data *) 
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Main Program 



begin 

page; 

initialize ; 

reset(data/ data ell*); 
reset(datam, * data magi*); 
reset(datac, *data curl*); 

rewrite(doutl, *data outl*); (* data after combine process *) 
rewrite(dout2, * tdata outl*); (* data of boundary *) 

for c : = 1 to col do 
for r : = row downto 1 do 

readln( data ,points( . r,c. ). el); 
for c : = 1 to col do 
for r : = row downto 1 do 
begin 

r eadln( dat am , t empmag) ; 
if tempmag < low_bound then 

points(. r ,c. )• i^^gtype := level 
else if tempmag > low_bound then 
points( . r , c. ) ‘ r^i^gtype := safe 
else points(. r ,c. )• r^^^gtype := unsafe; 
points( . r , c. ) . mag := tempmag; 

end; 

for c : = 1 to col do 
for r : = row downto 1 do 

readln(datac,points( . r ,c. ). cur); 



.V) 



(*printmagcur*) (* print magnitude and curvature *) 
rowlink; (* link points in row and make subregions *) 

printsub(points); (* print subregions created by rowlink *) 
collink; (* link points in row and column and make subregions *) 

calc_coeffs; (* calculate plane expressions for the 

subregions created *) 
combine( points , subregion) ; 

(* combine adjacent regions if they 
have similar plane expressions *) 
calc_coeffs; (* calculate plane expressions for the 

final subregions *) 
get_boundaries( subregion) ; 

(* get boundaries of subregions in linked list *) 
printsub2( points ) ; (* print subregions created by combine *) 

write_planar_data( points ) ; (* write out the elevation of the model *) 

write_boundary_data( points); (* write out the elevation of 

the boundaries of the subregions *) 



end. 
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APPENDIX C. 



PASCAL CODE FOR SMOOTHING ELEV ATION DATA 



AND CALCULATING GRADIENT 



(^Ss300000'0 

program calc_slope_curv_smooth( input, output); 

Author = Yee, Seung Hee. 

Date = 3 June 1988. 

Input = 1. Elevation data acquired from Fort. Hunter Liggett 
data base (40 by 40 square). 

Output = 1. File of magnitude of gradient of the points. 

2. File of curvature of the points (optional). 

3. File of smoothed elevation data (optional). 

Computer = Waterloo Pascal on IBM 370/3033AP, 

Naval Postgraduate School. 

Description = 

This program contains two major procedures: 
calc_slope_curvature and smoothing. 

The procedure calc_slope_curvature calculates the slope 

and curvature of each points, the procedure smoothing processes the 

elevation data to get smoothed elevation data. 

Original elevation data is represented in feet and there are 40*40 
data cells in one window. The distance between each cell isl2.5 
meters , so it needs to be converted to get a metric system. 



const 

distance_unit = 1; (* 3.28084*12.5 = 41.0105 *) 

row = 3; 
col = 3; 

curv_thresh = 0. 05; (* threshold for eigenvalues of biquadratic *) 
upbound = 0.6; (* slope > 60 % then unsafe slope *) 

lobound = 0. 1; (* slope < 10 % then level slope *) 

type 

elrange = 0. . 10000; 
pointrec = record 

el : elrange ; 
slope: real; 
cur: char; 
end; 

point_array=array ( . 1. . row, 1. . col. ) of pointrec ; 

var 

points ,spoints : point_array ; 

dat a , outs lope , out cur , smout el , smouts lope , smoutcur : text; 

counter , r , c , i , num: integer; 

k2,k3 ,k4 ,k5 ,k6 , slope, laml , lam2: real; 



procedure init(var pt: point_array); 
var 
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r,c,i,j: integer; 

begin 

for r : = 1 to row do 
begin 

pt( . r, 1. )• slope := 999; 
pt( . r,col. )• slope := 999; 
pt( . r , 1. )• cur : = ; 

pt( . r , col. )• cur := *u*; 
end; 

for c : = 1 to col do 
begin 

pt( . 1 ,c. ). slope := 999; 
pt( . row, c. ). slope ;= 999; 
pt(. l,c. ). cur : = ’u*; 
pt( . row, c. ) . cur := *u*; 
end; 



end; 

procedure calces lope_curvature(var pt: point_array); 
var 

r,c,i,j: integer; 

temp 1 , temp2 , t ermx , t ermy , dtermx , dtermy , 
templam , termxx,termyy,termxy,term : real; 
begin 

for r: = 2 to row-1 do 
for c: = 2 to col-1 do 
begin 

term : = 0.0; 
termx : = 0.0; 
termy : = 0.0; 
termxx : = 0.0; 
termyy := 0.0; 
termxy := 0.0; 

for i :=-l to 1 do 
for j : =-l to 1 do 
begin 

term := term + pt( . r+i , c+j . ) . el; 
termx := termx + j*pt( . r+i , c+j . ) . el; 
termy := termy + i*pt( . r+i,c+j. ). el; 
termxx := termxx + j*j*pt( . r+i , c+j . ) . el; 
termyy := termyy + i*i*pt( . r+i , c+j . ). el; 
termxy := termxy + i*j*pt( . r+i , c+j . ). el; 
end; 

k2 := termx/( 6*distance_unit); 

k3 := termy/( 6*distance_unit); 

k4 := (termxx/2 - term/3)/distance_unit; 

k5 := termxy/(4*distance_unit) ; 

k6 := (termyy/2 - term/3)/distance_unit; 

laml := ((2*k4+2*k6) + 

sqrt( abs( sqr( -2*k4-2*k6) -4*( 2*k4*2*k6-sqr(k5) ) ) ) )/2; 
lam2 := ((2*k4+2*k6) - 

sqrt( abs( sqr( -2*k4-2*k6) -4*( 2*k4*2*k6-sqr(k5) ) ) ) )/2; 
if abs(laml) < abs(lam2) then 
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begin 

templam := laml; 
laml := lam2; 
lara2 : = templam; 

end; (* eigenvalues of bigger slopenitude becomes laml *) 
slope := sqrt(sqr(k2) + sqr(k3)); 

pt(. r,c. )• slope := slope; 

if ( laml>curv_thresh) and ( lam2>curv_thresh) then 
pt(. r,c. )• cur := 'd* (* depression *) 

else if ( laml<-curv_thresh) and ( lam2<-curv_thresh) then 
pt( . r 5 c. ) . cur : = * h* (* hill or peak *) 

else if ( abs( laml)<=curv_thresh) and 
( abs( lam2)<=curv_thresh) then 
pt(. r ,c. )• cur := ^p* planar *) 

else if ( laml<-curv_thresh) and ( lam2>curv_thresh) then 
pt( . r ,c. )• cur := 's* (* saddle 

else if ( laml<-curv_thresh) and 

(abs( lam2)<=curv_thresh) then 
pt(. r ,c. )• cur := *r* (* ridge *) 

else if ( abs( lam2)<=curv_thresh) and 
( laml>curv_thresh) then 
pt( . r , c. ) • cur := 'v* (* valley *) 

else if ( lam2<-curv_thresh) and ( laml>curv_thresh) then 
pt(. r ,c. )• cur := *a* (* pass *) 

end; 

end; 



procedure smoothingC var points ,spoints: point_array); 
var 

^jC,i, j ,temp,dif : integer; 

begin 

for r : = 2 to (row-1) do 
for c : = 2 to (col-1) do 
begin 

if (points(. r,c. )• slope > lobound) 
then begin 

temp := 0; 

for i := r-1 to r+1 do 

temp := points( . i,c. )• cl + temp; 
for j := c-1 to c+1 do 

temp : = points(. r, j. )• cl + temp; 
spoints(. r ,c. )• cl := trunc( temp/6); 

end 

else spoints(. r ,c. )• cl := points(. r ,c. )• cl ; 
end; 

for r : = 1 to row do 
begin 

spoints(. r,l. )• cl : = point s(. r,l. ). el; 
spoints(. r ,col. )• cl := points(. r,col. )• cl; 
end; 

for c : = 1 to col do 
begin 

spoints( . 1 , c. ) * cl ; = points ( . 1 , c. ) . el; 
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spoints(. row,c. ). el : = points (. row,c. ). el 
end; 

end; (* procedure smoothing *) 



procedure check(var pt: point_array) ; 
var 

count ,countb,countp,counth,countd : integer; 

counte,countr, countv, count a : integer; 
r,c,aa,cc ; integer; 
begin 

count : =0; 
countb : =0; 
countp : =0; 
counth : =0; 
countd : =0; 
counte : =0; 
count r : =0; 
countv : =0; 
counta : =0; 
aa : = 0; 
cc : = 0; 

for r : = 1 to row do 
for c : = 1 to col do 
if pt(. r,c. )• slope < lobound then 
aa : = aa + 1 

else if (pt( . r,c. )• slope = 999.0) then 
countb : = countb + 1 

else if (pt(. r,c. ). slope > upbound) then 
cc : = cc + 1 



else 

case 



^t( . r ,c. ). cur of 



p 


countp 


1 = 


countp 


+1 


d’ 


countd 


1 = 


countd 


+1 


h^ 


counth 


: = 


counth 


+1 


e ’ 


counte 


: = 


counte 


+1 


r ’ 


countr 


; = 


countr 


+1 


v^ 


countv 


; = 


countv 


+1 


a’ 


counta 


: = 


counta 


+1 



end; 

writeln( 

writeln( 

writeln( 

writeln( 

writeln( 

writeln( 

writelnC 

writeln( 

writeln( 

writelnC 

writelnC 

end; 



level points are^ ^aa); 
unsafe points are ,cc); 
for safe points : ' 
count planar =' , countp); 
count hill =' , counth); 
count depression , countd); 
count etc =’ , counte); 
count saddle =', counts); 
count ridge =^,countr); 
count valley , countv); 

count pass =^, counta)*) 
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(*: 



Main Program 






begin 

reset(data, * data ell*); 
rewrite(outs lope, * data ma^l*); 
rewrite(outcur , * data curl^); 
rewriteC smoutel , * smdata ell*); 
rewrite(smout slope, ’ smdata raa^l* ); 
rewrite( smoutcur , * smdata curl^); 

for c : = 1 to col do 
for r : = row downto 1 do 

readln(data,points( . r,c. ). el); 
init( points); 
init(spoints); 

calc_s lope_curvature( points ) ; 

(* check( points) *) 
for c : = 1 to col do 
for r : = row downto 1 do 
begin 

writeln(outslope,points( . r,c. ). slope); 
writeln(outcur ,points( . r,c. ). cur); 
end; 

for r: = 1 to 10 do begin 

smoothingC points , spoint s ) ; 
calc_s lope_curvature( spoints ) ; 
check(spoints); 
points := spoints; 
end; 

for c := 1 to col do 
for r : = row downto 1 do 
begin 

writeln( smoutel , spoints( . r,c. ). el); 
writeln(smoutslope,spoints( . r,c. ). slope); 
writeln(smoutcur,spoints(. r,c. ). cur); 
end; 

end. 
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